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Abstract 

Four-fermi models in dimensionality 2 < d < 4 exhibit an ultra-violet stable renormalization group 
fixed point at a strong value of the coupling constant where chiral symmetry is spontaneously bro- 
ken. The resulting field theory describes relativistic fermions interacting non-trivially via exchange 
of scalar bound states. We calculate the 0(1/Nf) corrections to this picture, where Nf is the 
number of fcrmion species, for a variety of models and confirm their renormalizability to this order. 
A connection between renormalizability and the hyperscaling relations between the theory's critical 
exponents is made explicit. We present results of extensive numerical simulations of the simplest 
model for d = 3, performed using the hybrid Monte Carlo algorithm on lattice sizes ranging from 8 3 
to 24 3 . For Nf — 12 species of massless fermions we confirm the existence of a second order phase 
transition where chiral symmetry is spontaneously broken. Using both direct measurement and 
finite size scaling arguments we estimate the critical exponents /?, 7, v and 5. We also investigate 
symmetry restoration at non-zero temperature, and the scalar two-point correlation function in the 
vicinity of the bulk transition. All our results are in excellent agreement with analytic predictions, 
and support the contention that the 1/Nf expansion is accurate for this class of models. 
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1. Introduction 

The existence of a continuum limit of four-fermi theories is a long standing problem that goes back to 
the work of Nambu and Jona-Lasinio [1]. The original motivation for studying these models was the fact 
that they give a qualitatively good description of dynamical chiral symmetry breaking in strong interaction 
physics. The only degrees of freedom involved are fermions which interact via a short range interaction. 
As such four-fermi theories do not have a renormalizable perturbation expansion in powers of the coupling 
constant. However, at strong couplings, scalar bound states appear in the spectrum as a consequence of 
chiral symmetry breaking and the fermions interact by exchanging these composite scalars. The appearance 
of the scalars and pseudoscalars is tied to the Nambu- Goldstone realisation of chiral symmetry. By symmetry 
arguments, it can be seen that the emerging low-energy theory is similar to the linear u-model, with scalars 
coupled to fermions via a Yukawa interaction, which is known to be perturbatively renormalizable. In four 
dimensions, the equivalence between the two models was noted very early on [2] and has been exploited by 
many authors [3,4]. Meanwhile, it was discovered by Gross and Neveu that in two dimensions the model 
is asymptotically free and has a nontrivial continuum limit near the origin [5]. Wilson has argued that 
four-fermi theories are nontrivial for all d < 4 [6]. His arguments are based on the leading order results of 
the 1 /Nf expansion, where Nf is the number of fermion species in the model. Many other studies have been 
done using this expansion technique [7]; a series of articles by Guralnik et al. [3] indicated that below four 
dimensions, an equivalence between four-fermi and sigma models could be established. 

For definiteness, let us introduce the model we shall be dealing with in the bulk of this article. It is often 
called the Gross- Neveu model in the literature, and is the simplest relativistic theory of interacting fermions. 
The continuum space-time Lagrangian (we work in Euclidian space throughout) reads, 



where the index j labels the Nf species, each of which is described by a four component spinor. For purposes 
of analysis as well as computation, it is useful to introduce an auxiliary field a so that C becomes quadratic 
in tpji 



Gaussian functional integration over the a field in (1.2) restores the original Lagrangian (1.1) - however, it 
also serves as an interpolating field for the scalar bound state governing the fermion interactions described 
above. In the chiral limit fermion bare mass m — > the vacuum expectation value of a becomes a dynamical 
fermion mass and is a convenient order parameter for any chiral symmetry transition in the theory. Note 




(1.1) 




(1.2) 
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that in this limit Eq. (1.2) has a discrete Z 2 chiral symmetry: a — a; tpj 75^; an d i>j l— * ~"0j75- It 
is this symmetry which is spontaneously broken at strong coupling. The universal critical indices of this 
critical point have been calculated to 0(1 /Nf): for d = 3, they read 

8 8 18 16 

in the conventional notation of classical statistical mechanics [8] . These indices are particularly interesting 
because they are far from free-field Gaussian model indices. Apparently the critical field theory describes 
a strongly interacting fermionic system, complete with composite states, an interesting S-matrix, etc. The 
critical indices Eq. (1.3) lie at the heart of this article. The deviation of the critical indices from the Gaussian 
model indices are related to the anomalous dimensions of the various composite fields in the model, as we 
shall now discuss. 

Below four dimensions the Yukawa theory is superrenormalizablc with an ultraviolet fixed point at the 
origin and an infra red fixed point at strong couplings, where we recover scale invariance with scaling governed 
by non-vanishing anomalous dimensions. At this point the scalar degrees of freedom become composite and 
the low energy theory is the same as that of the four-fermi theory in the UV region. The compositeness 
condition is manifested in the vanishing of the scalar kinetic term in the continuum limit. Whereas at weak 
coupling (d^a) 2 is marginal irrespective of d, in the vicinity of the infrared fixed point the field a scales 
with an anomalous dimension which, because of unitarity, must be positive. This makes the kinetic term 
irrelevant at the IR fixed point. Since the IR limit of the Yukawa model now has the same relevant operator 
content as the UV limit of the four-fermi model, the latter must be renormalizable. Recently, several groups 
have studied the three dimensional model and demonstrated that Wilson's ideas survive beyond the leading 
order in 1/Nf. In particular, a series of papers by Rosenstein et al [9] have pointed out that the 1/Nf 
expansion around the fixed point in d = 3 is renormalizable. This result was extended to general d G (2, 4) 
in [10]. Explicit 0(1 /Nf) calculations appear in references [8,11]. Ref.[12] contains a rigorous proof of 
renormalizability of the three-dimensional theory. 

We should mention the physical meaning of the compositeness condition Z = 0, where Z is the wavefunc- 
tion renormalization factor associated with the scalar kinetic term, and its relationship with the fixed-point 
condition. Consider the general case of an interacting theory where a bound state \B > with binding energy 
Eb = —B < appears. The wavefunction renormalization constant is defined by 

Z = J2\<b\B>\ 2 (1.4) 

b 

where \b > stands for the bare (elementary particle) states. Using standard techniques, it can be shown that 
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Z satisfies the following equation [13] 

with G 2 (e) being the total decay rate of the state \B > proportional to some effective coupling constant. 
From Eq. (1.4), we see that Z = means that the bound state has no projection on the space of bare states. 
This is the compositeness condition. From (1.5) it is transparent that this condition puts an upper limit on 
the effective coupling G(e) thus yielding fixed-point condition. 

The renormalizability of a given theory can be argued in several ways. One is the old-fashioned way of 
constructing a finite theory by adding counterterms and ensuring that all divergences cancel order by order 
in some kind of expansion. Another is to start with a discretized, and therefore UV finite, theory and look 
for the region in the bare parameter space where a macroscopic length scale appears. The two methods are 
equivalent and we will use them both in this paper. The general idea of renormalizability is that the cutoff 
dependence can be absorbed into a finite set of bare parameters in such a way that the low energy physics 
which emerges is insensitive to the cutoff. Once this is done, it is possible to find lines of constant physics 
- renormalization group (RG) trajectories - in the bare parameter space. These lines are uniquely defined 
no matter which observable is taken. This way, the low-energy quantities depend on each other and not on 
the cutoff - we will see that this is possible if certain consistency conditions known as hyperscaling relations, 
governing the scaling of different physical quantities, are obeyed. 

The hyperscaling hypothesis [14] claims that the only relevant scale in the critical region is the macroscopic 
correlation length £. If this hypothesis holds, it is possible to do dimensional analysis using this correlation 
length as a scale. As a consequence, all dimensionless quantities should be independent of £. In the context 
of a simple model of ferromagnetism, the hypothesis can be stated as 

F sing = t 2 - a F(h/t A ) cx C d (1-6) 

where F sing is the singular piece of the free energy density and a, A = (36 are specific heat and gap exponents 
respectively. The external symmetry breaking field is h and t = T — T c is the deviation from the critical 
temperature. All thermodynamic quantities can be obtained by taking the derivatives of the free energy. 
In particular, the order parameter (magnetization), defined as < 4> >= dF sing /dh, satisfies the equation of 
state (EOS) 

< >= t 2 - a - A F'(h/t A ) = t^F'(h/t A ) (1.7) 

with the magnetic exponent defined as (3 = 2 — a — A. Similarly, the susceptibility exponent, 7, is defined 
via 

X = = t?- A F"{h/t A ) = t-<F"{h/t A ) (1.8) 



i.e. 7 = A — [3. The behavior of the correlation length in the scaling region is given by 

£ = t-»g(h/t A ) (1.9) 

Using the hyperscaling hypothesis (1.6) we relate aio v: dv = 2 — a. The exponents (3 and v can be related 
to the scaling dimension of the field (p. If a field <f>, which develops nonvanishing vacuum expectation value 
< <j> >^ 0, has scaling dimension = ^(d — 2 + 77), then hyperscaling implies < cf> >^ £ *. The vanishing 
of the order parameter is given by < (j> >^ ^ from where it follows that (3/v — d^. The other scaling 
relations are derived similarly. For a given channel, the corresponding correlation length (= inverse mass) 
may be defined by 

"S- I > (L10) 

Since hyperscaling is an important statement, we should explain its meaning and outline possible implica- 
tions. It is generally believed that the violation of hyperscaling leads to triviality. This is due to inequalitites 
between certain combinations of critical indices [15]. In the simplest cases like scalar field theories and spin 
systems, the quantity that measures the violation of hyperscaling is the dimensionless renormalized coupling. 
It is defined through the nonlinear susceptibility 

y (ni) / n d 3 < 6 > 

9 * = -??> = (L11) 

The renormalized coupling is essentially the properly normalized connected four-point function. In terms 
of the correlation functions, the nonlinear susceptibility is the zero-momentum projection of the four-point 
function 

X inl) = I < (Kx)</>{y)<f>(z)<f>(0) > CO nn (1.12) 
J xyz 

The normalization factors can be understood by recalling that susceptibility and mass are related to the 
wave function normalization constant Z through Z = x/£, 2 ■ I n the numerator we get one power of Z 1 / 2 for 
each field, which is thus compensated by % 2 : the factor of £ d accounts for the fact that there is one extra 
spatial integration in the numerator of (1.11). Both an d X can b e obtained by differentiating the free 
energy (1.6): x = t 2 - a - 2A F"{y) and x ( ™° = t 2 - a - 4A F lv {y). Together with the EOS for the correlation 
length (1.9), they give the expression for the renormalized coupling 

9R = t dv -^H(h/t A ), (1.13a) 

or, using the definition of the exponents (3 and 7, Eqs. (1.7,8), the alternative form 

g R = t- 2A+ - 1+d "H(h/t A ) (1.136) 
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In specific (ferromagnetic) models an estimate of the exponent can be made. Since multi-spin correlations 
cannot extend over a larger range than pair correlations, the following inequality holds [15]: 



2A < 7 + dv 



(1.14) 



We can trade t with the correlation length in Eq.(1.9) and rewrite Eq.(1.13) in terms of £ as 



9r = ^ 2A -^- dv)/v H(h/t A ) 



(1.15) 



Since gn is dimensionless, hyperscaling implies that it must be independent of £ and the exponent must 
vanish. In this case the renormalized coupling is a function of only one bare variable. Strict inequality in 
Eq.(1.14) implies triviality. An equivalent way of stating the above inequality is dv > 2 — a. It implies that 
the singular part of the free energy vanishes no faster than the prediction of hyperscaling i.e. F S i ng > £~ d . 

In a similar fashion the scaling of mass ratios can be derived [16]. If hyperscaling is satisfied, for any 
particular pair of masses that obey Eq.(1.9) e.g. M v , M a , we have 



where G(y) is a universal function. Comparing it with the expression for the renormalized coupling, we see 
that both observables depend on the same variable. One of the relations, R = G(y) or gn = H(y) can be 
inverted to solve for the bare variable e.g. y = H^ 1 (gn). This defines an RG trajectory for each value of the 
renormalized coupling: h = H^ 1 (gjj)t A . This is then used to obtain the relation between the two observables 
R = R{gii)- The same manipulation can be done with two mass ratios. We note that the important point 
in the inversion is that both observables depend on just one bare variable, so that the inverse relation can 
always be found, at least in some regions of parameter space. This would be difficult to achieve otherwise, 
and is a consequence of hyperscaling. 

To summarize, the requirement of having a macroscopic correlation length as the only relevant scale leads 
to relations between the critical exponents. These relations are 



All the other relations are obtained from Eqs.(1.17) using the definitions of the critical exponents e.g. 
2[35 - 7 = dv, 2/3 + 7 = dv. 

Next we review briefly how these ideas are realized in four-fermi theories, at least to leading order in 
l/Nf. We shall see in the next section, where detailed calculations are presented, that sublcading corrections 




G{h/t A ) 



(1.16) 




(1.17) 
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do not change the physics qualitatively. The critical exponents for the model (1.1) have been calculated in 
Rcf. [8]. The physical fermion mass M satisfies a gap equation, 

M = m- g 2 <ipip >, (1.18) 

In the chiral limit M is proportional to the order parameter < tpip >. In the large Nf limit, only the simple 
fermion tadpole contributes to < ipip >. The exponents (3,8 can be obtained from the gap equation. It is 
convenient to use the definition of the critical coupling 1 = 4g 2 J^l/q 2 , which is cutoff dependent, to rewrite 
the gap equation in the form 

r j\f3 

tM + m = 4g 2 —-- —r- (1.19) 

where t = (g 2 — g 2 )/g 2 . Below four dimensions, the integral is IR divergent in the absence of the fermion 
mass M. It can be evaluated exactly and it is of the order of M d ~ x . Thus, in the limit when either bare 
mass m or t vanish, we have 

(1.20) 

M-m 11 ^-^ (t = 0). 
Since < >~ the exponents are: (3 = l/(d— 2), 5 = d— 1. 

The exponents v, 7, t] are obtained from the scalar propagator D a . To leading order in 1/Nf it is given 

by 

D- 1 (k 2 ) = l + g 2 Il(k 2 ), (1.21) 

where II(fc 2 ) = — tr S(q + k)S(q), and S is the fermi propagator (ifi + M)^ 1 . We associate the scalar 
mass with the inverse correlation length 1/M 2 = £ 2 . The susceptibility and wavefunction renormalization 
constant for the scalar field are: x 1 = 1 + .9 2 n(0), Z~ x = g 2 n'(0), which using (1.10) leads to 



_J dD-\k 2 ) 

D-\k 2 ) dk 2 



k 2 =0 



The scalar mass and renormalized coupling defined in this way are thus 

M 2 - /y- 1 - 1+ g 2n (°) ~ 2 - Za 2- 1 M 2V 

where the renormalized Yukawa coupling gn is dimensionful . It has dimension (mass)^ 4_d ^/ 2 and, since it 
is a low-energy quantity, it depends on the physical mass M a . The dimensionless coupling is thus 

The central question in the study of four-fermi theories is whether the bound states remain composite or 
become pointlike in the continuum limit. An important consequence of compositeness is the non-triviality 
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of the theory; ie. the dimcnsionlcss interaction strength remains non-vanishing in the continuum limit. The 
physical reason in the case of the linear cr-model is easy to understand. The generic coupling between the 
fermions and scalars is given by the Goldberger-Treiman relation: g v fj — M/f^, where M is the fermion 
dynamical mass and f„ is the pion decay constant related to the pion radius by rv ~ 1/ f v - This way, the 
Yukawa coupling becomes 9V// ~ Mr n and vanishes in the limit of pointlikc pions. In other words, since the 
scalar mass is not protected from receiving large corrections (of the order of the cutoff scale) via radiative 
corrections, the effective interactions will always be zero-range in the continuum limit. Such a theory does 
not have a physical scale and must be trivial above two dimensions. The only way to have a nontrivial limit 
is if the exchanged particles are composite. We will see later that this requirement is achieved when the bare 
parameters are tuned so as to make all power-law divergences disappear. 

We restrict the discussion to 2 < d < 4 in what follows to avoid violations of scaling although, regarding 
the flow of RG trajectories, the conclusions will be the same for d = 2 (but not for d = 4). Using the gap 
equation, we get (the quantity Il(fc 2 ) will be evaluated in detail in the next section) 

i m 12 o -, 9 m 12 9 , bq 2 

* =M + —) MZ > M ^M Z+ —) M > Z =M^' 

where b = 2T(2 — d/2)(d — l)/3(4n) d ^ 2 . We measure all the masses in units of the momentum cut-off. Note 
that the wavefunction renormalization constant Z vanishes in the continuum limit M — > 0: this is precisely 
the compositeness condition [13]. In the chiral limit m — > the susceptibility is x _1 oc M d ~ 2 , which implies, 
using Eq.(1.20), that oc t, i.e. 7 = 1. Similarly, the a mass in Eq.(1.25) scales as M a oc M oc t d ~ 2 giving 
v = l/(d— 2). The exponent rj is extracted from the power law decay of the scalar correlation function at 
the critical point. Simple power counting in Eq.(1.21) gives 

which gives 77 = 4 — d. 

To obtain the RG trajectories and relate the low energy observables, we use Eq.(1.25) to calculate the 
mass ratio R = M 2 /M 2 and the renormalized coupling: 

1 m 12 2 _ 1 

R - bg^M^ + Jd~Ty 9r ~ bR 2 -m- (1 - 27) 

It is clear from the second relation that the lines of constant mass ratio are also lines of constant renormalized 
coupling. Regarding the ratio between the scalar and fermion masses we observe two things. In the limit 
d — > 4, since 1/b oc 4 — d, R — > 4, which is the expected result if the scalar is a weakly bound fermion 
anti-fermion state. For arbitrary d in the chiral limit, we obtain the result R = R c = 12/ (d — 1). At first 
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sight this seems to contradict the interpretation of a as a weakly bound state - however, this result is due 
to our using the definition (1.22) for the scalar mass, together with the fact that the numerical coefficient of 
Z is far from unity. As we shall confirm in the next section, the scalar propagator has a pole at k 2 = —AM 2 
in the chiral limit for all d. Away from the chiral limit, however, D a no longer has a pole, which means that 
the definition (1.22) is more useful from our perspective. Secondly, consider the gap equation; 

cg 2 M d - 2 = t+^ (1.28) 

with c = 8r(2 - d/2)/(d - 2)(4 7 r) <i / 2 . To leading order in 1/N f this can be regarded as the EOS. From the 
expression for R and the gap equation we eliminate M to obtain 

M-fAj (129) 

where A — /35 = (d — l)/(d — 2). This is the expected scaling of the RG trajectories m ~ t A . The lines of 
constant R fall into two subfamilies. For t > 0, R — R c < c/b and t < 0, R — R c > c/b. In the chiral limit 
R = R c . As the binding increases, the mass of the composite, measured in units of the constituent mass, 
decreases. 

We comment about the theory above four dimensions where hyperscaling is violated and where, as a 
consequence, the continuum limit is trivial. It is easy to see that, for d > 4 the gap equation (1.19) is regular 
in the M — > limit. As a consequence, the corresponding exponents /?, S are the same as in four dimensions; 
(3 = 1/2, S = 3. Again, the reason for this is the IR behavior of the theory. The gap equation in this case is 

Cg 2 A d - 4 M 3 = tM + m, (1.30) 

where we display the explicit cutoff dependence, since the physical mass will turn out not to be the only 
scale. It is easy to verify that the expression (1.29) for the mass ratio will remain unchanged. However, the 
renormalized coupling will change; 

d-4 

B V A 



g 2 R = ZM d -Y = I ( ^f) . (1.31) 



The wavefunction renormalization constant is obtained as the k 2 coefficient in the scalar propagator. To 
leading order it is 

which is clearly finite in the M — > limit for d > 4. The relations for the remaining critical exponents 
follow from Eqs.(1.25): 7 = 1, v = 1/2, r\ = 0. For d > 4 hyperscaling is clearly violated: indeed with 



mean field values of the exponents Baker's inequality (1.14) becomes d > 4. Two things are clear: firstly 
the trajectories of constant R and gn do not coincide; and secondly in the limit M/A — > 0, the renormalized 
coupling vanishes. Because of the violation of hyperscaling the theory is trivial. 

In four dimensions, hyperscaling violations are logarithmic and are not reflected through the violation of 
the relations between the critical exponents. Explicit calculation to leading order gives 



So, the theory is trivial, but the critical exponents satisfy hyperscaling relations. The reason for this is that 
the scalars are pointlike: Z~ x is logarithmically divergent in the UV cutoff. 

In section 2 of this article we present a calculation of 1/Nf corrections to the model at length, for 
arbitrary dimensionality 2 < d < 4. We will explicitly demonstrate the renormalizability of the 1/Nf 
expansion to next-to-leading order in this regime, and also calculate the critical indices and show that 
they continue to satisfy the hyperscaling relations. We will argue that the physical assumptions underlying 
renormalizability and hyperscaling are equivalent, and demonstrate this using the chiral symmetry of (1.1,2). 
For completeness we also generalize our results to models where the spontaneously broken symmetry is 
U(l) <g> U(l) or SU(2) ® SU(2). Some of these results were presented in briefer form in Ref. [8]. 

In the rest of the article we shall fix d = 3 and study the model by computer simulation methods. There 
are several motivations for doing this. First, we can investigate the validity of the 1/Nf expansion by a truly 
non-perturbative numerical scheme. We shall see that standard methods of extracting critical behavior from 
computer simulations yield results in good agreement with the 1/Nf expansion when Nf is chosen large - 
Nf will be set to 12 in our numerical work. Second, we can study the model for small Nf, where it may 
have some relevance to the strongly correlated electron dynamics underlying high T c superconductivity, and 
see if there are dramatic changes as Nf decreases. Although we have some results at small Nf, this article 
will concentrate on the first objective and extract a number of results from Nf = 12 simulations which 
provide support for the usefulness and reliability of the 1/Nf expansion. Our main results are a numerical 
verification of the critical index predictions of Eq.(1.3) to leading order. 

After describing the lattice formulation of the model and the hybrid Monte Carlo algorithm in Sec. 3, we 
turn to numerical results. Most of our simulations are performed directly in the chiral limit m = 0. In Sec. 
4 we consider the model for Nf — 12 on symmetric lattices ranging in size from 8 3 to 20 3 . Measurements of 
the vacuum expectation value So =< & > clearly show a chiral transition in the lattice model at an inverse 
coupling 1/g 2 near unity. Measurement of T, and its susceptibility \ as a function of the coupling 1/g 2 allow 




Sir 2 



(1.33) 
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us to measure the critical indices (3 and 7 directly. The leading order terms in Eq. (1.3), /? = 7 = 1, will be 
confirmed with good precision. The size dependence of the critical point will also provide a rough estimate 
of u, using finite size scaling arguments, but more quantitative results will be presented in later sections of 
this article. In particular, the model will be considered at nonzero temperature in Sec. 5 by simulating it on 
asymmetric lattices, N T x TV 2 with N = 36 lattice spacings and N T ranging from 2 through 12. We shall find 
a second order phase transition at a critical temperature in good agreement with large 1/Nf calculations. A 
measurement of the critical index v will follow. In Sec. 6 we introduce a small bare fermion mass m into the 
theory's action and measure the critical index 8 at the zero temperature critical point as well as at the finite 
temperature transition. The first set of measurements, done on 8 3 , 16 3 and 24 3 lattices, are nicely consistent 
with the leading order term in Eq.(1.3), S = 2. The second set of measurements uses a 12 x 36 2 lattice 
and gives 8 T = 3.0, the expected Gaussian model index for the temperature-driven phase transition [17]. In 
Sec. 7 we turn to measurements of the two-point scalar correlation function, and our attempts to extract 
the scalar mass using the formalism developed in Sec. 2. As we shall see, this requires some novel fitting 
techniques, and the results presented here are probably only a first step towards a complete understanding. 
Our results are summarized in Sec. 8. 



2. Calculations of 0(1/Nf) Corrections 

In this section we will calculate the lowest non-trivial corrections, ie. those of 0(1/Nf), to the model in 
the vicinity of the strongly-coupled fixed point at g = g c . We will work in dimensionality d where 2 < d < 4, 
and have chosen to follow the treatment of Gat et al [11], and define the Lagrangian involving the auxiliary 
scalar field in the chiral limit as follows: 



where a sum on i = 1, . . . , Nf is understood. Notice that in this section of the paper the a field is normalized 
differently to that of Eq.(1.2); since a is auxiliary this has no physical consequences. The constants Z$, 
Z /y , and g are cutoff-dependent, and must be adjusted at each order of the 1/Nf expansion to keep Green 
functions finite, the relevant ones being the fermi propagator Sf, the scalar-fermion vertex T^^, and the 
scalar propagator D a . It should be noted that the constants Z^, Z a are chosen dimensionless for convenience 
and are not directly related to the constant Z of the previous section, which has dimension (mass) 2 . For 
additional simplicity we work for the most part in the chiral limit m = 0. In the broken phase t > 
we will choose renormalization conditions so as to keep the physical fermion mass, that is, the pole of the 




(2.1) 
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renormalized fermion propagator, a finite constant independent of 1/Nf. 
Our Feynman rules are thus: 

4 0) - 1 =Z v ,(^+S ); 
r (0) _ _ 9 



f 



(2.2) 



where the dynamical fermion mass S appearing in the propagator is defined by the expectation value of 
the scalar field (this is equivalent to the definition of S used elsewhere in the paper): 

^ gzl 



/N f 



< a > . 



(2.3) 



We determine Eo from the equation of motion 8C/5a = 0, yielding the gap equation. To leading order, 



9 2 N f 



which immediately leads to an expression for g in terms of Sq and the cutoff A: 



1 



9 2 (4rr)*(d-2) 



A d ~ 2 

rXf) 



s^- 2 r(2-f) 



(2.4) 
(2.5) 

of expression (2.5) 



Details of how J is defined for d € (2,4) are given in an appendix. The limit Sc 
defines the critical coupling g c : for d = 3 wc obtain g c — 2A/ir 2 . To the same leading order the scalar 
propagator D a is now given by the sum over all numbers of chained fermion - antifcrmion bubbles n(fc 2 ): 
D-\k 2 ) = Z a -Il(k 2 ) 
= Z a g 2 tv 



1 



1 



(2.6) 



= ^.9 ; 



2 5-rfr(2- |) 



(fc 2 + 4S 



Jo 



dxx 2(1 -a;) 2 2 ; 



(47r)5 Jo 

where use has been made of relation (2.4). The integral over x in (2.6) defines an incomplete Beta function, 
which in turn may be expressed in terms of a hypergeometric function. The result is 



° {k) ~ a9 ~^ji sp - F(1,2 2,21 m y 



For d = 3 we recover the usual form [9] 



D-\k 2 ) = Z a g 



2 (fc 2 + 4S 2 ) tan-^V^^Eo) 



(2.7) 



(2.8) 



27rVfc 2 

Notice that for fc 2 <C Z? CT resembles the canonical form for a boson field: (fc 2 + 4£q) ; indeed, in the 
chiral limit it has poles at fc 2 = — 4Sq. However, for fc 2 » Eg, the relevant limit for the calculation of 
0(1 /Nf) loop corrections, the asymptotic form of (2.7) is quite different: 

1 f A d \ 



lim D a {k 2 ) 



k 2 — *-oo 



z*g 2 \(k 2 )l-\ 
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(2.9a) 



where 

A d = 3 S— s • (2.96) 

To renormalize the model at leading order, we may take Z$ = 1, and the physical mass M to be equal 
to So- At higher orders this will not be the case. Finally to render D^ 1 finite we define 

The exact value of the right hand side is unimportant provided it is finite, since it cancels from all calculations 
of physical quantities. In effect we have specified a to have dimension [|], but since it is an auxiliary field 
this is purely a matter of convention. 

The simplest corrections at 0(1/Nf) are one loop corrections to Sf and T a ^, as shown in Fig. 1 (in 
fact, there is a two loop contribution to at this order [9]: since, however, it involves a closed fermion 

loop with an odd number of legs, it is finite and can be ignored). The fermion self-energy T,(p) is given by 

m - E - x ± x ^ + ; ) + So ZW). (2.11) 

The integral over q is logarithmically divergent but straightforward, and leads to the result for the full inverse 
propagator 



. (d-2) C d A \ / C d A 

^( 1 + ^TAV ln Mj +So {'-m-^M 



(2.12) 



where we have used the fact that M — S to leading order, and the numerical constant Cd, which we will 
find to be ubiquitous, is given by 

Cd= T(2-l)T(l)B(l,l-iy (2 ' 13) 
In d = 3, Cd = 4/-7T 2 . For finiteness we thus require 



and 



V - ! Ml-^^ln^). (2.1.-,) 



Similarly, at vanishing external momentum we find the vertex to 0(1/ Nf): 

9 r, i A Cd , A 



= —*=Zi,ZZ [ 1 - ^ In - ) . (2.16) 



/]V/ V 2N f M 

The requirement that this expression be finite gives a condition on the combination Z a g 2 : 

( 1+ 2§^) ; 
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le. 



Z a g 2 = 



1 



1 , 2(d-l) C dl A 

1 + ^^]v7 ln M 



(2.17) 



As we saw in Eq.(2.7), there is a second condition on Z a g 2 at leading order from the requirement that the 
scalar propagator be finite. This continues to be the case at higher orders, as we shall see: consistency of 
the two conditions at each order is vital for the renormalizability of the model. 

To fully specify the renormalization constants to 0(1/Nf) we must calculate g as a function of physical 
mass M using the two-loop gap equation (Fig. 2). The 0{1/Nf) contribution is: 



s 1 



,(0)2 



4 

'W f 



(3p 2 -2p. g -Sg) 2 



4(d- l)A d - 



2 d-i^- 2 r(i-|)(rf-i) 



(2.18) 



JV / (47r)*r(|)(d-2) 
using expression (2.9) for _D CT . The full expression for g is thus: 



(4rr)i 



8A d - 



(47r)«r(f)(d-2) 



1 - 



2iV> 



8S 



ci-2 



(47r)»(d-2) 



r(2-f) 1 



>d-3 



(d-1) 



7V f 



(2.19) 



Once again, we extract the critical coupling in the limit E — > 0: 



8A d ' 



(47r)ir(|)(d-2) 



1 - 



2A^ f 



(2.20) 



We have now used up all our freedom in adjusting the constants g, Z^p and Z a by absorbing divergences 
in the fermion self-energy and vertex corrections, and by relating quantities to a physical scale via the gap 
equation. Power counting shows, however, that we have not yet exhausted the list of divergent graphs at 
0(1/Nf): there remain the two-loop contributions to the scalar propagator shown in Fig. 3. If the model is 
renormalizable, then the 0(1/Nf) expresssion for Z^g 2 already derived in (2.17) together with the definition 
of g via (2.19) must suffice to cancel the UV divergences in D a . We will show this by explicit calculation in 
the symmetric phase, where the massless nature of the fermion makes the loop integrals much simpler. Of 
course, the UV divergence structure of a field theory should be insensitive to the particular phase in which 
it resides. We begin by following the renormalization convention of [10]: the leading order expression for the 
inverse scalar propagator in the symmetric phase is now (Cf. (2.6)) 



D- 1 (k 2 ) = Z <7 g 2 
The integral over q is easily done, to yield 

D-\k 2 )=Z a g 2 



g 2 



8A fl 



1 



1 



(2.21) 



g 2 (An)Hd-2)T(i) 



(2.22) 
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where Ad is given in (2.9). We now impose the off-shell renormalization condition 

££V) = ^. (2-23) 

where /i is some arbitrary scale parameter. The power-law divergence in (2.22) can now be removed, first 
by specifying 

Za9 2 = f- d -, (2.24) 

then, substituting (2.22) into (2.23) we find 

1 8A d - 2 u d - 2 

(2.25) 



9 2 (An)Hd-2)T(l) A d 
With Z a and g 2 both fixed, D^ 1 is completely finite: 

D ° l{k2) = A^( {k2)i ~ 1 + ^ 2 )- (2 - 26) 

We see that in the symmetric phase D a {k 2 ) has no poles on the physical sheet, and so \x cannot be interpreted 
as a scalar mass: indeed the a field in this case does not interpolate a stable particle. Instead, [i can be 
regarded as the width of an unstable resonance, but as shown later, it can still serve as an inverse correlation 
length in the statistical mechanics sense. The critical coupling g c is obtained in the limit \i — > 0, whereupon 
we find the same result as for the So —> limit of (2.5). Note that in (2.25) g < g c , which is consistent with 
being in the symmetric phase. To leading order no further divergences are present, so as before Z$ = 1. 

It is important to note that the asymptotic form of (2.26) matches that of the propagator in the broken 
phase (2.9), as it must do. It is therefore no surprise that we obtain almost identical results for the fermion 
wavefunction and vertex corrections at 0{1/Nj) as for the broken phase (Fig. 1), the only difference being 
that now no mass renormalization is required, due to the chiral symmetry of the model. We obtain two 
conditions for the renormalization constants essentially identical to Eqs.(2.14,17): 

(d — 2) Crf A 

,2_ 1 f-t I 2(ri - 1) C'd . A 
fi d - 2 V d N f ii 

with Cd given by (2.13). 



= 7^ 1 + ^-T^ Tf:^-l (2-28) 



Since there are no non-trivial solutions to the gap equation for g < g c , the renormalization of the model 
can only be completed by extracting g(A, /j,) from the full inverse scalar propagator to 0(l/Nf), which 
necessitates a two-loop calculation. The relevant contributions are shown diagramatically in Fig. 3 (once 
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again, there is a finite three-loop contribution which we ignore [9]). For external momentum k, the two- loop 
contribution to D~ 1 (k 2 ) is 



A, 



q ( 9 2 )f- 1 +^-2 



1 



1 



+ 2tr 



1111 

ij> i{j> + i)ij> i{j> + $) 



(2.29) 

We will spare the reader the full details of the calculation, and merely remark that the integral over q is 
considerably simplified by our choice of momentum routing through the internal scalar line. The divergences 
are generically of two kinds: a power-law form 

A d-2 



(a — 2) [i 



(2.30) 



and a momentum-dependent logarithm 



(ii) ~ (fc 2 ^- 1 ^ 



A 



{ M |A}' 



(2.31) 



where the symbol {fi\h} denotes a function of both \x and k which simplifies in three limits of interest: 



limjtt | k} = v A: 2 ; 
lirn{/x | fc} = /x; 
| im 2 {M I fc} = A 4 - 

Further details on how forms (i) and (n) are extracted from (2.29) are given in the appendix. 
After some work we arrive at the full expression for D~ 1 (k 2 ) to 0(1/Nf): 



(2.32) 



D~\k 2 ) = Z a g 2 



+ 



8A 



d-2 



(47r)S(d-2)r(j 



8S(|-l,|-l)r(2-|) ^ d 



(4^)^(1) 



N 



f 



— — rr-M 2 In- --(fc) 2 In-— — ^ 
(d-2) fi) d V 



= Z a g 2 



^_^_ ( ,_ 1) ^ ln A + (^f 1 _^i)^ ln a >^ 



(2.33) 

where g c is identical to the result for the broken phase, Eq.(2.20). Once again, we renormalize using condition 
(2.23), plus the expression (2.28) for Z a g 2 coming from vertex renormalization. We find an expression for g: 



d-2 



g 2 g 2 c + A d 



i + 



(d-\){d-2)C d . A 



d 



, r ln- 

N f m 



(2.34) 
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Substituting back into (2.33), we find the renormalized inverse scalar propagator: 

D-Hk^-Z aA^~" (l 2(rf-l) Q ln A\ (fc 2 )^ 1 ( 2{d-l) C d A \ 



(2.35) 



1 



V d N f n ) 



A df i d 

This expression is manifestly UV finite both for fc 2 = 0, and in the critical limit [i — > 0, although in the 
latter case there is the usual IR divergence associated with massless particles. 

Since we have been able to adjust Z a and g to eliminate all divergences in Green functions to 0(1/JV/), 
at least in the symmetric phase, we have succeeded in our aim of showing the model is renormalizable to 
this order. In effect the fact that the combination Z a g 2 serves both to correct the fermion-scalar vertex 
and the inverse scalar propagator is really a statement that the renormalized four-point fermion Green 
function is finite to this order (Fig. 4). If we recall that the underlying model we started from was one 
with an explicit four-fermi interaction, we can now understand why Z a g 2 is apparently overdetermined: in 
the four-fermi model there can only be two adjustable parameters: the unphysical Z^; and the strength 
of the four-point coupling. Another consequence is that for k 2 ^> M 2 ,/i 2 , logarithmic contributions to the 
four-point scattering of the form ln(fc/M) must also cancel, which means that in this asymptotic limit the 
renormalized four- fermion scattering amplitude assumes a universal form A d /N fk d ~ 2 . This receives no 1/iVj 
corrections if the model is renormalizable. 

We now consider the model from the viewpoint of statistical mechanics, by calculating the critical ex- 
ponents associated with the chiral symmetry breaking transition, and demonstrating a relation between 
renormalizability and hyperscaling. All essential information about a continuous phase transition is encoded 
in its critical exponents, and it is straightforward to calculate these to 0(1/Nf) from the gap equation. The 
order parameter is of course the condensate < -0V >j but from (2.4) we can equally well consider So, since it 
shares the same non-analytic behavior as g — > g c +, A/M — > oo. Parameterizing the distance from criticality 
by t = g 2 (g^ 2 — 9~ 2 ), we hnd from (2.19) the critical exponent [3: 

l 

t oc oc S^ 2 , (2.36) 

whence 

"=(^2) +0 (^)' (2 - 37 > 

To determine the other critical exponents describing the critical behavior of the order parameter, we 
need to consider the effects of introducing an explicit bare mass term Z^m'tpitpi into the Lagrangian (2.1). 
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To leading order M = So + m , and we obtain an inverse scalar propagator which now has a a contribution 
diverging as a power of A at leading order: 



D-\k 2 -m) = Z a g\r 



i 



i 



i 



2Y{2 - |) (k 2 + AM 2 ) 



(47T)i 



M 4 - 



F(l,2 

A d-2 



{4n)i(d-2) S 



2 ' 2 ' 4M 2 



(2.38a) 



— - M d - 2 T(2 - | ) 



fc 2 



lim D-\k z ;m) = Z a g 2 



(k 2 



r(f) 

m / A d ~ 2 



- Af d - 2 r(2 - 1) 



(2.386) 



A d (47r)i(d-2)£ V r (f) 
It is interesting to note that away from the chiral limit D a no longer has poles in the complex fc-plane - 

therefore to define a "cr mass" we are forced to use the definition (1.10,22). If we substitute the asymptotic 

form (2.38b) into the two-loop gap equation (2.18), with M replacing S in the fermi propagators, we obtain 

the full gap equation in the presence of external masss: 



g 2 9 2 



8M d - 2 T(2- |) 



8(d- 1) 



C d A d - 2 m, (. (d-2)C d Z 



(4rr)*(d-2) (4tt)t (d - 2) 2 r(|) N f S 
Setting g — g c and taking the limit ra«Eo~Mwe find 



In 1 + 



m 



m 
772 



(d-1) C d So 
(d - 2) TV/ n m 



8Sg _1 r(2- 
(47r)*(d-2) ' 



le. 



(d-1) 



1 + ^ 

AT/ 



(2.39) 



(2.40) 



(2.41) 



where we have assumed that the logarithm exponentiates beyond 0(1/Nf). To extract 7, we must return to 
the full expression (2.38a) for D a (k 2 ;m), substitute it into (2.18) and evaluate d/dm\ m=0 on the resulting 
gap equation. After some work we find 



Cd, A\ 



7V + {d - 1) w f l ^ )=^ + ^ 

Near criticality \ 3> 1, from which we can infer 



^d-2 



8r(2 - 1) 

(4tt)7 



using (2.36), ie: 



X oc [ d - 2 +( d - 1 ) c 'd/ Ar /] 0, t -[l+(d-l)C d /(d-2)N f ] 



(d-l)Q 

7 (d-2)A/ 



(2.42) 



(2.43) 



(2.44) 



We now see that our derived values for (3, S and 7 are consistent to 0(1/Nf) with the scaling relation 



1 = 0(6-1), 
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(2.45) 



which raises the issue of whether the exponents calculated in this way obey scaling or hyperscaling relations. 
To develop this picture we need to calculate further critical exponents not directly determined by the behavior 
of the order parameter, and hence need to supply some physical insight. The obvious choice for correlation 
length is the inverse physical fermion mass M _1 . Using Eqs. (2.15) and (2.36) we find 

t oc C b oc mV-W-V-WW, (2.46) 

leading to a prediction for the exponent v. 

i r. fd-nci (2 47) 



(d-2) 



1 . l (d-l)C d 



d N f 



We can calculate the final exponent rj in two different ways. Firstly and more directly, consider the 
renormalized scalar propagator (ie. two-point correlator) (2.35): at criticality we expect the scaling form 

D- 1 oc fc 2 -". (2.48) 

From the [i — > limit of (2.35) we deduce 

Of course, strictly we should consider correlations of the scalar composite iptp, which defines the true order 
parameter. Following the treatment of [10], we can renormalize this local operator, and use (2.4) and (2.15), 
to write 



M = - "' J 1 * < V^i >= ~ * J 1 " Z U < {MiU >, (2.50) 



N f N f 

where the constant has been introduced to renormalize V>,Vi- Now, both sides of (2.50) are cutoff- 
independent (since M is a physical mass), as is (V^VOft itself; note also that a factor of Z^, is absorbed in 
taking the trace to calculate < tpip >. We conclude that 

A |v ( Z mV ^) = 0. ( 2 - 51 ) 



ie. 



- _ A d _g 2 Ad (Z M \_ g 2 (d-l) C d 

where we have used (2.5,15), and the quantity denoted by 7^ is the anomalous dimension of the operator 
-iptp defined in the conventional field-theoretic sense [18]. At the fixed point g = g c , we can use the resulting 
■y^ — (d— 2) + (d— l)Cd/dNf to derive r\. We expect the renormalized two-point correlation function 
< (Tpip(0))n{\pTp{x))R > to scale asymptotically in momentum space as 

T^ 2) (fc 2 ) oc fc- [d - 2 4> 2 ^ oc fe-P-^**], (2.53) 
18 



where we have introduced the canonical (ie. free-field) scaling dimension for tptp: 

d^ = d-l- (2.54) 

or in other words rj = d—2j^ using (2.48). The value for 77 thus derived is identical to (2.49), demonstrating 
that these rather formal arguments about operator renormalization are valid, and that for all practical 
purposes a and ipip can be considered identical. 

The calculation of the critical exponents is now complete: to 0(l/Nf) the values we have computed for 
/?, S, 7, v and 77 satisfy the hyperscaling relations (1.17): 

2/3 + 7 = dv; 2f3S - 7 = dv; (3 = \>{d - 2 + rj); 7 = 1/(2-77). (2.55) 

In the realistic case d = 3 discussed in the rest of this paper, the values coincide with those given in Eq.(1.3). 
It is also interesting to examine the limits d — > 4_ and d — > 2+. Noting that Cd vanishes in each case, as 
(4 — d) and as (d — 2) respectively, we find for d — > 4 we recover mean-field exponents; 

P = b 5 = 3 ; 7 = 1; ^ = 5; ( 2 - 56 ) 

while for d — ► 2 we find 

It is interesting to note that we could have inferred the rcnormalizability of the model merely by con- 
sidering the exponents, without performing the integral (2.29). Suppose we had assumed hyperscaling to 
derive 77 from the gap equation exponents (3, 5, and 7, and made no further assumptions about operator 
renormalization, then we could have used relations (2.48) to predict the divergence structure of the scalar 
propagator at criticality: 



2{d-l) C d A 
d N f n k 



(2.58) 



The argument of the logarithm is determined on dimensional grounds, since the cutoff is the only other 
scale left in the problem. Similarly we could use an alternative definition of susceptibility, together with the 
derived values of 7 and u, to deduce the momentum-independent divergence: 



X" 1 = D~ 1 (k 2 = 0) oc f oc - M d - 2 



x 2(d-l) C d]n A 



(2.59) 



d N f M _ 

In both cases we can check that Z a g 2 given in (2.17) eliminates all dependence on the cutoff A and renders 
D^ 1 finite. Hence we can infer renormalizability at this order without ever doing a two-loop calculation 
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(although the gap equation has two loops, it has one fewer fermi propagator, and no external momentum 
dependence, and so is much simpler!). 

There would appear to be an intimate connection between hyperscaling and renormalizability which 
we can demonstrate to all orders in 1/Nf without doing any detailed calculations. Consider the effect of 
shifting the a field in the Lagrangian (2.1) so that the fermions pick up a mass term Z^EoV'iV'i- It is then 
straightforward to derive the following functional identity [19]: 

d S 2 T v/777 s 3 r 



<9S Si^Stp ggi SipSipSa' ' 



(2.60) 



where a' is the shifted scalar field and V = T[ip, ip, a] is the effective action. We can use this relation in the 
broken phase to relate the fcrmion self-energy T,(k 2 ) to the full fermion-scalar vertex [9]: 

E(0) = -Eo-Vpi-I^O), (2.61) 
gZ3Z^ 

The chiral symmetry of the original Lagrangian implies E(0) is proportional to E , and power counting 
shows that and hence E are at most logarithmically divergent functions of A. Using the definitions 

implicit in (2.11-16), we deduce 

Eo 1 

Z iP Z M s o oc - — Z a 2 , 

ie. 

Z 2 M cx Z a g 2 . (2.62) 
The constant of proportionality is 0{M 2 ~ d ) from (2.10). Now, we know that Zm has the form 

Z M -l + aln^, (2.63) 

i i 

where we assume a ~ 0(1 /Nf) is small, so that So oc M l ~ a . By definition, we also have t oc Eg oc M" , 
which enables us to relate a to (3 and v assuming that higher powers of the logarithm exponentiate as required 
by the hypothesis of power-law scaling near a critical point: 

l-a=^, (2.64) 



so from (2.63,64) we deduce 



Z-,,r x | + 2(1- ^) In A (2.65) 



From previous we know that if the model is renormalizable Z„g 2 must suffice to render the inverse scalar 
propagator D^ 1 finite, by cancelling both fc-dependent and fc-independent logarithmic divergences, once 
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power-law divergences have been removed by tuning g 2 to its critical value (Cf. Eq.(2.33)). Consider each 
case separately. At criticality, we must have on dimensional grounds 



D-\k 2 ) oc fc 2 -" ~ k a 



1 - (4 - d - 77) In ^ 

K 



(2.66) 



If the logarithmic divergence in (2.66) is to be cancelled by the Z a g 2 of (2.65), then 

2(l-§) =4-d-rj, 

ie. 

/3 =ii/(d-2 + »j), (2.67) 
which is one of the hyperscaling relations (2.55). Similarly in the limit k 2 — > we have 

A" 



so that renormalizability implies 



(2.68) 



2/3 + 7 = di/, (2.69) 

which is another hyperscaling relation. Hence with a few simple assumptions about the form of divergences 
at 0(l/Nf) and beyond, we can see the equivalence of hyperscaling and renormalizability. This puts the 
renormalizability of the model on a very physical footing; indeed, the main hypothesis behind both properties 
is the assumption that there is only one important physical length scale (correlation length) whose divergence 
in cutoff units at the critical point controls the scaling of all other physical quantities. In this model the 
physical length scale is the renormalizcd fermion mass M in the broken phase and the scalar width \x in the 
symmetric phase: both ratios A/M and A/fj, diverge with the same critical exponent v in the chiral limit at 
the phase transition, which defines the continuum limit. The continuity of v is demanded on physical grounds 
[20] , since otherwise the free energy of the system would be discontinuous across the phase transition. 

Suppose we had found that in the broken phase the value of v demanded by hyperscaling differed from 
the value vm derived by assuming that the renormalized fermion mass M is indeed an inverse correlation 
length. We would then have been forced to conclude that cither the interaction strength tended to zero as 
.9 ~~ * 9c {vm > v), so that the model described free massive fcrmions, or that the fermions became very 
massive and decoupled [i>m < v), leaving a theory of non-interacting scalar bound states. In either case the 
model would have had a trivial continuum limit. However, this scenario is impossible in the present model, 
since we know from the identity (2.60,61) that the definition of renormalized fermion mass is tied to the 
renormalization constant Z a g 2 , which in turn cancels out divergences in the inverse scalar propagator if and 

21 



only if hyperscaling is obeyed. Therefore a model in which um 7^ v would be non-renormalizable. Thus in 
this model, and as we shall see, in similar four-fermi models, the issues of hyperscaling, renormalizability, 
and non-triviality are inseparable. 

Finally, we note that in addition to references [9,11], there have been other calculations beyond leading 
order in the Gross-Neveu model for 2 < d < 4. Our result for Z$ is consistent with the calculation of 
Hikami and Muta [21], modulo different definitions of the trace over the Dirac algebra. We also agree with 
the results of [9,11,22] for d = 3. Gracey [23] has also computed the exponent r\ (in our notation) and has 
actually obtained Z$ to 0(1 /N 2 ): his sophisticated methods rely on being exactly at the critical point and 
hence cannot be used to calculate the other exponents. Other recent calculations beyond leading order for 
general d have appeared in [24]. 

Finally in this section, for completeness, we consider extensions of the Gross-Neveu model to cases where 
the spontaneously broken symmetry is a continuous one [5]. There are two interesting cases; the symmetries 
in the Lagrangian are either U(l) £ ® U(1) JJ (= U(l) y ® U(1) A ): 

C = $iMi--^ [(V^i) 2 - (^75^) 2 ] ; (2.70a) 

or SU(2) L ® SU(2) fl : 

g 2 - 

C = ipipfiipip - [(ipip^ip) - (ip ip l 5 T pq ^ lq ) 2 ] , (2.706) 

where f are the Pauli matrices, normalised to tr(r a r' 3 ) = 2<5 Q/3 , and we show the indices p, q running from 

1 to 2 explicitly In 2 < d < 4 dimensions we define 75 thus: 7! = 1, {75, 7^} = 0, tr(7 5 7 Ml . . . 7^ n ) = 0. We 

can immediately bosonize the models, and introduce renormalization constants as in (2.1): 

C = Z^ipitfipt H j=ZipZ^ [atpiipi + i7r^75^i] + ^-Z^cr 2 + 7r 2 ); (2.71a) 

V N f ' 1 

C = Z^iWi + -l= z 4> z l [aWi + in-AlbTipi] + ^Z^a 2 + t? 2 ), (2.716) 
where a is an auxiliary scalar and 7r an auxiliary pseudoscalar. Note that in each case the combination 
4> = a + in (a + iir) is proportional to an clement of the chiral group, so that eg. in the SU(2)®SU(2) 
case we have ipL l_ * UipL; ipu 1— > VipR; <fi 1— > V^U -1 , where ipL,R = (1± 75)^/2, and U and V are SU(2) 
matrices. This property does not in general extend to higher flavor groups such as U(n). Model (2.71b) 
closely resembles the original Gell-Mann - Levy a- model [25]. Notice also that the same renormalization 
constant Z^ serves for both a and 7r fields: this is the result of a Ward identity which we discuss below. 

To leading order in l/Nf model (2.70a) has an identical gap equation to (2.5) and expression for D a 
(2.7). The only new feature at this order is the pion propagator, which in the broken phase is given by 

D^(k 2 ) = ^ 2 ^-#-^F(l, 2 - i ■ >; -^). (2.72) 
(47r) 2 L 
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The pole at k 2 — explicitly shows that the ir field has a Goldstone mode in the broken phase. In the 
symmetric phase we find as before (2.24,26): 



D-\e) = D-\e) = ^[{e)i-^+^]. 



(2.73) 



The SU(2)<g)SU(2) model goes through in very similar fashion apart from an overall factor of 2 in the gap 
equation from the trace over Pauli indices, and a corresponding factor of \ in D a and D^, the latter of 
course now carrying an implicit 5 af3 . 

The calculation of 0(1/Nf) corrections proceeds as before, the only new feature being the appearance of 
a logarithmic divergence in the 0(1/ Nf) correction to the gap equation containing an internal tt line. The 
renormalization constants are found to be, for U(l)<g)U(l): 



Z^f, — 1 — 



Z M = 1 + 
Z^g 2 oc 1 + 



(d - 2) C d 



d 


N f 


(d- 




d 


N f 


2(d- 


-2) C d 



A _ 
A 



(2.74) 



d N f 



In 



M ' 



with g 2 given by the gap equation: 



8A fl 



(47r)*r(f)(d-2) 



(d-2) 
N f 



8Sg- 2 r(2 - f) 
(47r)*(d-2) 



2 d ~ 2 (d-2) 
Nf 



(2.75) 



Similarly for SU(2)®SU(2) we find: 



7 -1 ( d - 2 )^, n A - 
iV7 ln M' 

= 1 H In — 

2d N f M 

7 2 1 , (d-4) g d A 



(2.76) 



and the gap equation: 



1 



16A d - 



(47r)ir(f)(d-2) 



1 



(2d - 5) 
2N f 



16S^ 2 r(2- I) 



2 d - 2 (2d-5) 3^ A 
2JV, + 2 (d ~ 2) AV lD ^ 



(2.77) 



(4tt)S (d-2) 

Once again, we find that the combination Z^g 2 derived via vertex renormalization is also sufficient to cancel 
divergences in the two-loop scalar and pseudoscalar propagator corrections. The calculation of critical indices 
is straightforward, and the results for these, plus the critical coupling g 2 , are tabulated in Table I. It is easily 
checked that the exponents satisfy the hyperscaling relations (2.55) to 0(1/Nf). 
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The arguments (2.60 - 69) showing the relation between renormalizability and hyperscaling go through 
exactly as before, except that now because the broken symmetry is continuous, we can start from a fully- 
fledged Ward identity, eg. for U(l)(g>U(l) we have 

- <<j>T^{ P T T = 0)= l -{S F \ lb }. (2.78) 

Using relations (2.3) and (2.12) we arrive at a result analoguous to (2.62): 

Z 2 M cx Z^g 1 . (2.79) 

To complete the argument we need a further Ward identity: 

< a > r ff(T ^(p w = o) = |{r CT ^, 75 }. (2.80) 

The four-point function Y is constructed from superficially convergent graphs, so needs no renormal- 
ization. Assuming 75 commutes with the scalar vertex, we see that Z a = Z n = Z^ as promised. Equivalent 
arguments hold for the SU(2)<g>SU(2) model. 

We end this section with a comment about proving renormalizability. We have stressed the physical 
equivalence of renormalizability and hyperscaling for model field theories of this kind, and, we suspect, for 
a variety of other models with a non-trivial fixed point in 2 < d < 4 [6,26]. A crucial point in our argument 
has been the existence of a chiral symmetry leading to identities such as (2.62,79). The model will remain 
renormalizable so long as this symmetry is only broken by soft terms, ie. a bare fermion mass. (The effects of 
introducing a symmetry-breaking term (-tpiipi) 3 into C and the consequent non-renormalizability of the model 
have been discussed in [22]). It has been stated [9,22] that the renormalizability of the model follows from 
power-counting considerations alone, once the leading order form of the scalar propagator (2.7) has been 
established. We feel that this ignores the necessity for the consistency relations between the renormalization 
constants we have presented here - hence the chiral symmetry of the model must play a fundamental role 
in the full proof. Of course, this is not a new situation: power-counting alone does not suffice to prove the 
renormalizability of gauge theories. Moreover, the SU(2)<g>SU(2) model provides a specific example where 
next-to-leading order corrections render scalar exchange harder than leading order, since 77 > 4 — d. Power- 
counting arguments to analyze the degree of divergence of the graphs would fail in this case unless the 
momentum dependence of the dressed vertices is taken into account. The demonstration that 0(1 /Nf) and 
higher corrections to scaling dimensions cancel in the correct fashion, so as to leave the set of primitively 
divergent graphs given by the leading order predictions unchanged, depends on consistency relations derived 
from Ward identities. 
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3. Lattice Formulation of the Gross-Neveu Model 

The Gross-Neveu model in its bosonized form (1.2) may be formulated in three dimensions on a spacetime 
lattice using the following action: 

Nf/2 ( 1 \ AT 

s = E E^(^.»xi(i/) + gE^(*)xiW E + d>I> 2 (*)> 

i=l \ a;,y a; <£,ie> / y £ 

where Xi?Xi are Grassmann- valued staggered fermion fields defined on the lattice sites, the auxiliary scalar 
field a is defined on the dual lattice sites, and the symbol < x,x > denotes the set of 8 dual lattice sites x 
surrounding the direct lattice site x. The lattice spacing a has been set to one for convenience. The fermion 
kinetic operator M is given by: 

■M x ,v = \ E l S v,x+ti - S y,x-(i\ , (3-2) 

where rf^(x) are the Kawamoto-Smit phases (— l) Xl ^ 

The Gross-Neveu model in two dimensions was first formulated using auxiliary fields on the dual sites in 
reference [27]. We can motivate this particular scheme by considering a unitary transformation to fields u 
and d [28]: 

A (3.3) 

Here Y denotes a site on a lattice of twice the spacing of the original, and A is a lattice vector with entries 
either or 1, which ranges over the corners of the elementary cube associated with Y, so that each site on 
the original lattice corresponds to a unique choice of A and Y. The 2x2 matrices Ta and B A are defined 

by 

(3.4) 

Ba = (-t 1 ) a H-t 2 ) M (-T3) A3 , 
where the r M are the Pauli matrices. Now, if we write 

iron = (3.5) 

and interpret q as a 4-spinor with two flavors counted by the latin index a, then the fermion kinetic term of 
the action (3.1) may be recast in Fourier space as follows: 

d 3 k 



Skm = / (2^3 EE ! j&WCV ® h)qi(k) sin2fc M + %(fc)(74 <g> r*)^(fc)(l - cos2fe M )|, 



(3.6) 



I (J, 
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where 

{i,U= ( T " T ) ; (74)a/ J =f. 1 ; (7 5 W=f 1 l2 ) , (3.7) 

the second set of (2 x 2) matrices in the direct product act on the flavor indices, and the momentum integral 
extends over the range G (— tt/2, tv/2]. At non-zero temperature the lattice has finite extent in the 
temporal direction, and J dko is replaced by a sum over N T /2 modes, where N T is the number of lattice 
spacings in the time direction, and antiperiodic boundary conditions are imposed on the fermion fields. In 
the classical continuum limit lattice spacing a — > 0, the flavor non-diagonal terms vanish as 0(a), and we 
recover the standard Euclidian form q~j$qj, where the flavor index j now runs from 1 to Nf. 

The interaction term between the fermions and the auxiliary requires a little care. We label the dual site 
(x + i, x + 1, x + i) by (A;Y) where x corresponds to (A;Y). In this particular labelling the a fields are 
not all equivalent. For instance, for A = (0, 0, 0) it is easy to show the interaction term transforms to 

<T(0;Y)qi(Y){U®l2)qi(Y). (3.8) 

For A = (1, 0, 0), however, it is more complicated: 



(7(1; y) 



qi(Y)(l 4 ® U)qi(Y) + ^A+(qi(Y)(U ® U + hm ® r*)q t (Y)) 



(3.9) 



where A+f(Y) = f(Y + £i) — f(Y) ~ 2ad fl f(Y). We see that there is an extra term containing non-covariant 
and flavor non-singlet interactions, which is formally 0(a). If we used a formulation in which the a fields lived 
on the direct lattice sites, then such non-covariant terms would contribute at O(a ) [27]. It is straightforward 
to extend this analysis to all A; we conclude that the interaction term may be written 



2 (E a(A; Y^j g j (y)(l 4 ® l 2 )q t {Y) + 0(a), 



(3.10) 



which is clearly of the same form as the equivalent term in (1.2). In principle the 0(a) non-covariant 
terms could be important once quantum loop corrections are computed, and should be analyzed further. 
In the two-dimensional lattice Gross-Neveu model their effect is discussed in [29], where it is argued that 
in a perturbative expansion, even the on-site auxiliary formulation (with non-covariances present at O(a )) 
yields the correct physics, that is with no important continuum symmetries violated. This analysis appears 
to be readily extended to the 1/Nf expansion in the three dimensional case. 

Thus we see that the lattice action (3.1) reproduces the bosoniscd Gross-Neveu model, at least in the 
classical continuum limit. Most importantly, (3.1) has a discrete Z 2 global invariance under 

Xi (x) i * (-ir+ x ^ Xi (x); xi{x) ~ -(-i) I1+a!a+Xs Xi(aO; » (3-Ho) 
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ie. 

qi(Y)»(j 5 ®l)qi(Y); q^Y) ^ -qi(Y)(j 5 ® 1); a(x) ^ -a(x). (3.116) 

It is this symmetry, corresponding to the continuum form ip i— > 75"0, i/> l— * ""075 > which is spontaneously 
broken at strong coupling, signalled by the appearance of a non- vanishing condensate < \X > or equivalently 
< g(l4 (8> l2)g >■ As discussed in Sec. 2, an equivalent order parameter is the scalar field expectation S . 
To leading order in I/AT/, we can compute the fermion tadpole explicitly to yield the lattice gap equation, 
relating S to 1/g 2 (Cf. Eq.(2.4)). Due to the simplicity of the loop integral, the contributions from the 
0(a) terms in (3.10) vanish, and we find 

S 1 _ r /2 d 3 k -^[sh^fc^fg. l 2 ) + (l-cos2fc M )( 74 ®r;)] + S (l 4 ® 1 2 ) 

— r = Trtr^F = / —- — rrtr - — 5 . (3.12) 

5 2 V A/2 (2tt)3 sin 2 + 

Using standard techniques, we arrive at a form suitable for numerical quadrature: 

1= I™ dae-^lhl&), (3.13) 
9 Jo z 

where Jo is the modified Bessel function. A plot of S vs. 1/g 2 is shown in Fig. 5 - we see that in this 

regularisation scheme, to leading order in 1/Nf, the bulk critical point 1/g 2 ~ 1.011a -1 . 

The action (3.1) was numerically simulated using the hybrid Monte Carlo algorithm [30], in which the 
Grassmann fields are replaced by real bosonic pseudofermion fields (f>(x) governed by the action 




(3.14) 



where 

M xyi j = M X y5ij + 5 xy 8ij\ <t{x). (3.15) 

Note M is strictly real, and that in this model we are able to work directly in the chiral limit bare fermion 
mass m = 0, since the matrix M has non- vanishing diagonal entries, and can always be inverted. Integration 
over 4> yields the functional measure y/Aet(M t M) = detM if the determinant of M is positive semi-definite. 
This condition is fulfilled if AT//2 is an even number. 

The hybrid Monte Carlo algorithm works by evolving the fields through a fictitious "simulation time" 
r using a Hamiltonian H = it 2 /2 + S[u], where ir(x) is a momentum field conjugate to a. The fields are 
sampled after evolving through a period of fixed or variable r known as a trajectory. Detailed balance is 
maintained by accepting or rejecting the entire trajectory according to a Metropolis step calculated using the 
probability weight exp(— AH), and ergodicity by regenerating the {ir} after each trajectory using gaussian 
random numbers. The Hamiltonian dynamics is simulated by finite difference equations parameterized by 
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some fundamental interval dr: in the limit dr — * we expect the dynamics to conserve energy, and hence the 
Metropolis acceptance rate to approach one. The art of hybrid Monte Carlo lies in tuning the parameters 
of the simulation so that dr can be made as large as possible, to reduce the amount of computer time 
per trajectory, whilst maintaining a reasonably high acceptance rate. One method, as originally discussed 
in [30], is to alter the couplings of the "guidance Hamiltonian" used to evolve the system through the 
trajectory; only the "acceptance Hamiltonian" used in the Metropolis step need be exact. The idea is that 
couplings, masses, etc. may receive 0{dr) corrections, so that the algorithm is in effect producing exact 
Hamiltonian dynamics at a different point in coupling space. In the work presented here we have found that 
the acceptance rate of the algorithm may be optimized by tuning the parameter Nf in the action (3.14), 
so that the guidance Hamiltonian uses Nf > Nf. Interestingly, because quantum fluctuations are 0(1/Nf), 
the amount of computer time needed to achieve comparable statistical accuracy is roughly independent of 
Nf : the only requirement that increases with Nf is the storage of the psuedofermion fields. 

The hybrid Monte Carlo algorithm proved to run sufficiently quickly and efficiently that unusually quan- 
titative results could be obtained for large enough Nf. Considerable details about the runs will be given 
below as different computer experiments are described, but a few words about the parameters used in the 
algorithm should be recorded here. In order to maintain a high acceptance rate so the algorithm would 
produce configurations that explored phase space rapidly, we tuned both Nf and dr. Typically as the lattice 
size was increased dr had to be taken smaller and Nf approached Nf. For example, on a 6 3 lattice when 
simulating the Nf — 12 theory, the choices Nf — 13.35 and dr — 0.25 gave acceptance rates greater than 
93% for all couplings of interest. To maintain this acceptance rate on a 12 3 lattice we picked Nf = 12.38 
and dr = 0.125, and on a 24 3 lattice Nf = 12.17 and dr = 0.050. Measurements were taken only after a 
reasonable time interval had passed - it is pointless to take measurements after each sweep because sequential 
measurements are highly correlated, but measurements should not be taken so infrequently that information 
is lost. For most of our runs we chose a trajectory length r equal to half a time unit - in other words dr 
multiplied by the number of sweeps between refreshment was chosen to be 0.50. Note that on the 24 3 lattice 
a trajectory rises to 10 sweeps and the runs are proportionally more compute intensive than on small lattices. 
Of course, we never assumed that sequential trajectories were statistically independent. Wc used standard 
binning methods to calculate the variances in our measurements for each observable. In a typical run of 
5,000 trajectories, say, a list of 5,000 measurements of the order parameter, etc. would be made and those 
lists would be analyzed to find the number of truly statistically independent measurements. Near the critical 
point we found the usual symptoms of critical slowing down - tens of trajectories were needed to decorrelate 
measurements of observables such as the order parameter. Variances could be calculated after each run and 
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the run extended if greater accuracy were needed. In our tables of results to be discussed below, we typically 
list the number of trajectories used and the statistical error bars for our raw data sets. Runs as long as 
10,000 trajectories were performed near the sytem's critical point to extract scaling laws and critical indices 
with good accuracy and confidence. Runs of this length are ten times as long as those currently practised 
in lattice QCD. We have such good statistics here because of the lower dimensionality of our model and 
the fact that each sweep has much fewer floating point operations. In addition, the number of sweeps of 
the conjugate gradient routines needed to effectively invert the lattice Dirac operator during each sweep is 
much smaller than in lattice QCD. Typically, only 10 - 100 conjugate gradient sweeps are needed for each 
inversion, while in lattice QCD ten times as many are needed in the present generation of simulations. The 
reason for the difference presumably lies in the fact that Eq.(l.l) has only a discrete rather than a continuous 
chiral symmetry. We very carefully monitored the convergence of the various conjugate gradient routines 
we used. There was the conjugate gradient used in the guidance Hamiltonian, that used in the acceptance 
Hamiltonian, and that used in measurements of the chiral condensate. The stopping residuals were typically 
chosen per lattice site to be 10~ 4 , 10~ 6 , and 10~ 4 respectively. Since conjugate gradient routines are only 
approximate, they could introduce an unwanted systematic error into the algorithm. Therefore, we carefully 
checked that our observables were insensitive to the size of the stopping residuals. Since the conjugate gra- 
dient algorithm converges monotonically at an exponential rate, it is relatively cheap (in computer time) to 
choose the stopping residual very safely and conservatively. Lengthy test runs were made to assure ourselves 
that all was in order. 

We have not chosen to explore different values of Nf systematically in this study, but rather to thoroughly 
explore the system's critical behavior with the choice Nf = 12. However, in pilot studies on small lattices 
we did comparative runs at Nf — 6, 12, and 24. The results from a 12 3 lattice for the expectation value 
S =< a > vs, 1/g 2 are plotted in Fig. 5. We see that chiral symmetry is indeed spontaneously broken, 
and for 1/g 2 between 0.5 and 0.8 the measurements are in fair agreement with the leading order prediction 
(3.13): the points lie systematically below the line, which we interpret as being due to 0(1/Nf) corrections 
as predicted by Eq.(2.19). In Table II we give values for S (iV/ = oo) given by (3.13) together with the 
measured normalized deviations iV/(Eo(oo) — Xo(iVf)). The numbers for Nf = 24 and 12 are consistent 
within errors, implying that the deviation is 0(1/Nf). For Nf — 6 the deviation is consistently larger, 
possibly because for this small value 0(1 /N 2 ) effects are becoming significant. A similar trend was found 
in a high-statistics study of the two-dimensional Gross-Neveu model [31]; here Nf was varied between 4 
and 120, and the 0(1/ Nf) correction calculated exactly using a lattice regularization. It was found that an 
0(1 /N 2 ) term was required to accurately fit data for Nf < 12. More work is needed in three dimensions 
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before reliable conclusions can be drawn; for one thing, at Nf = 6 the model does not have a positive definite 
determinant, and runs at adjacent values of Nf will be needed before the results of Table II can be taken 
seriously. 

Finally we note that for Nf — 6,12, chiral symmetry is restored for 1/g 2 > 0.9, and for Nf = 24 
it is restored for 1/g 2 > 1.0. Finite volume effects to be discussed in the next section make the precise 
determination of the critical coupling 1/g 2 difficult. In the rest of the paper we will restrict our attention 
to the case Nf = 12, and concentrate on high precision studies in the region 1/g 2 > 0.7, using a much finer 
grid of 1/g 2 values. This will enable a quantitative description of the critical scaling properties of the model. 



4. Symmetric Lattice Simulations and the Bulk Critical Point 

Our emphasis in this section is the discovery of the critical point predicted by the large- Nf expansion and 
a numerical calculation of its critical indices, j3 (magnetic), 7 (susceptibility) and v (correlation length). This 
program will be very successful for (3 and 7 while an accurate determination of v will only occur once the 
behavior of the model at non-zero fermion density has been investigated [32] . To find the critical point we can 
simply calculate S , the vacuum expectation value of the auxiliary field. Since So is proportional to the chiral 
condensate < ipip > and since its vacuum expectation value spontaneously breaks chiral symmetry, it is a 
particularly convenient order parameter. On small lattices, however, vacuum tunneling processes which take 
S to — S can obscure the critical point. The same problem affects computer simulation calculations of the 
magnetization in the Ising model, for example, and can be controlled by adding a small symmetry breaking 
field to the action. We chose not to do that here, although later we shall use this trick to measure the critical 
index 5. Instead, we monitored each computer simulation run for vacuum tunnelling events. Away from the 
critical point in the symmetry-broken phase such events were so rare that good measurements of S and 
its susceptibility \, given by the variance of So, were possible. There are presumably two reasons for this. 
First, we chose Nf = 12 which is sufficiently large that fluctuations and tunnelling processes were unlikely. 
And, second, our lattices were sufficiently large (8 3 - 20 3 ) that the probability of tunnelling events was 
highly suppressed. However, in the immediate vicinity of the critical point tunnelling events were sufficiently 
common that our computer data was not useful. Luckily we shall see that the scaling laws (critical indices) 
we seek were apparent in the data over regions of coupling where tunnelling was not a problem. 

In Figs. 6 - 9 we show the data for S vs. 1/g 2 and the reciprocal of the susceptibility l/\ vs. 1/g 2 
on both the ordered and disordered sides of the transition. Our raw computer data, statistical errors and 
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number of trajectories of the hybrid Monte Carlo algorithm used for each measurement can be found in 
Tables III - VI. For 1/g 2 < 1/g 2 , the critical point, we expect a non-analytic vanishing of S , 

(41, 

where (3 is the magnetic critical index. This simple power behavior is only expected for couplings 1/g 2 
sufficiently near the critical point 1/g 2 - In leading order of the large Nf expansion (3—1. Similarly the 
susceptibility should diverge at 1/g 2 ; 

X = A (\-\ y ) 9 2 > 9 2 c 

W c 9 2 J ^ (4 2) 

x=s (^~£i 7 gl>g2 

with the same index 7 above and below the transition. In leading order of the large Nf expansion 7=1. 

The data plotted in Figs. 6-9 are beautifully fitted by the large Nf predictions (3 — 1 and 7=1. 
We are unable to extract the small finite Nf correction (for Nf = 12, 7 = 1 + 8/N f n 2 + 0{Nj 2 ) ~ 1.068, 
(3=1 + 0{NJ 2 )) expected in this simulation. It may be necessary to work very close to the critical point 
to find evidence for 7 = 1.068 rather than 7 = 1, and that would require larger lattices to evade vacuum 
tunnelling apparent in the figures and the Tables III - VI which accompany them. Note that the error bars 
on the measured £0 values are smaller than the circles themselves of the figures. The error bars on the 
susceptibility \ are considerably larger (5 - 10%) because it measures the width of the fluctuations in the 
order parameter. In Fig. 10 we plot lnS vs. \n(l/g 2 20 — 1/g 2 ) where l/g 2 20 — 0.955 is the best estimate 
for the critical coupling of a 20 3 lattice inferred from Fig. 9. The straight-line fit in Fig. 10 gives (3 = 1.00 
and is in almost perfect agreement with the measurements at couplings 1/g 2 = 0.70 - 0.825. Comparing 
Tables V and VI we see that over these couplings the E measurements are in agreement to better than 1% 
on the 16 3 and 20 3 lattices. Closer to the critical point we do not have comparable control over finite size 
effects. It would be interesting to redo Fig. 10 on a larger lattice and see if the S values at 1/g 2 > 0.825 
approach the scaling curve. Similar curves for the susceptibility x can be made (lnx vs. \n(l/g 2 2Q — 1/g 2 )) 
and the critical index 7 found to be 1.0(1). Much greater statistics would be needed to reduce the error 
on the determination of 7 to one percent where the finite Nf corrections in Eq. (1.3) could be probed. A 
major barrier to our obtaining a higher precision determination of 7 is the uncertainty in the exact critical 
coupling and the fact that our estimates of it are hindered by vacuum tunnelling. 

One notices from Figs. 6-9 that the peak of the system's susceptibility shifts with the linear lattice size 
L. Finite size scaling arguments relate the size dependence to the correlation length exponent v; 

1 1 1 

-^ = aL-u, (4.3) 



9l(L) g 2 
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where l/g 2 (L) is the coupling where \ peaks on a L 3 lattice and l/g 2 is the infinite volume limit. Unfortu- 
nately the results shown in Figs. 6 - 9 do not determine l/g 2 (L) accurately enough to determine v well. In 
Table VII we give the l/g 2 {L) values with error bars and in Fig. 11 we plot Eq. (4.3) for v — 1 and infer 
l/g 2 — 1.00. In the next section of this article we shall study the four-fermi model on asymmetric lattices 
(nonzero temperature) and will achieve quantitative results for v. 



5. Asymmetric Lattice Simulations, the Critical Temperature, and the Index v 

We simulated the four-fermi model on asymmetric lattices, N T x N 2 , to determine its critical termperature 
T c where chiral symmetry is restored and to determine the order of this transition. Lengthy runs, several 
tens of thousands of trajectories of the hybrid Monte Carlo algorithm, in the vicinity of the phase transition 
on lattices with N T ranging from 2 to 12 and N set to 3N T showed no convincing evidence for metastability 
or tunnelling between a symmetric and an asymmetric vacuum. So, we concluded that the transition was 
second order in agreement with the large Nf analysis of this transition [17,32]. In Fig. 12 we show the 
histograms of measurements of £ (we reserve the subscript zero for the zero temperature value) in the 
ground state at couplings l/g 2 = 0.86, 0.865 and 0.870 on a 10 x 30 2 lattice. At 1/g 2 = 0.860 the system 
is in a chirally broken phase with £ very small, 0.039(1), but definitely nonzero. At l/g 2 — 0.865 vacuum 
tunnelling occurs between two states with £ = ±0.030, and finally, at l/g 2 = 0.870 the distribution of £ 
measurements is centered around the origin indicating symmetry restoration. For l/g 2 values decreasing 
below l/g 2 p c {N T = 10) = 0.865(5), the mean values of £ increase smoothly, which is indicative of a second 
order phase transition. Since T = l/N T = 0.10 in lattice units, we see that the critical temperature measured 
in units of £ is T c /£ = 0.64(4) in this case. (Note from Fig. 9 that £ = 0.157(1) at coupling l/g 2 = 0.865 
at zero temperature.) For values of l/g 2 greater than l/g 2 ^ c {N T = 10) the mean values of £ on the 10 x 30 2 
lattice were always consistent with zero and histograms of particular runs were well fitted with normal 
distributions centered at zero. Similar investigations on lattices with iV T ranging from 2 to 12 led to the 
estimates of the critical couplings l/g 2 jc (N r ) recorded in Table VIII. Using Fig. 9 to read off zero temperature 
vacuum expectation values of £ at these couplings produces the estimates of T c /Eo shown in Table IX and 
plotted in Fig. 13. 

These results should be compared to the large Nf prediction [17,32], 

| = 21^2 * 0J2 < 

which is plotted in Fig. 13. Our N T — 12 result is within a standard deviation of the exact result and the 
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trend of the simulation results, T c /E increasing slowly with N T , is quite satisfactory. 

Next we can use the data in Table VIII to obtain a better estimate of the correlation length index v. 
The shifts of the critical couplings l/g^ c (N T ) with lattice extent in the temporal direction should scale as 

where 1/g 2 is the large volume limit of the bulk critical coupling. From Fig. 10 we know that 1/g 2 is between 
0.95 and 1.00, but unfortunately cannot be specified with any greater precision. In Fig. 14 we plot Eq. (5.2) 
with the data of Table VIII for the choices of 1/g 2 — 0.950, 0.976 and 0.995. In each case linear fits of 
\n{l/g 2 — l/gp c (N r )) vs. lniV T are possible for N T = 6, 8, 10 and 12. The slopes determine the critical index 
v, and we have from the fits v = 0.81 if l/g 2 c = 0.950, v = 0.94 if 1/g 2 = 0.976 and v = 1.07 if 1/g 2 = 0.995. 
Clearly our precision for a good determination of v is limited by our uncertainty in the critical coupling. 
Nonetheless, our best estimate for v, 

v = 0.94(13), (5.3) 

is consistent with and close to the large N T prediction of unity. In the second paper of this series [32] where 
we study the four-fermi model at nonzero chemical potential a more accurate and independent determination 
of v will be given which lies within the error bars of Eq. (5.3). 



6. The Critical Index S for the Bulk and the Finite Temperature Transition 

The critical index S controls the response of a magnetic system's magnetization M at criticality to a 
small external magnetic field h: 

l 

M\ T = Tc ochs. (6.1) 

For the four-fermi model M maps into the chiral order parameter S ( or < V'V' >) an d h maps into a bare 
fermion mass (which explicitly breaks the chiral symmetry), so 

l 

Yi Q \ g 2 =g 2 ocm*. (6.2) 

For large Nf the index 5 is predicted to be 2, Eq.(1.3). We attempted to measure 5 by introducing a bare 
fermion mass into the theory's action, and simulating the system on 8 3 , 16 3 and 24 3 lattices at criticality 
but variable m. Since we do not know 1/g 2 exactly we did these simulations at 1/g 2 = 1.00 and 0.975. In 
addition, since we do not know beforehand how small m must be taken to see the scaling law Eq. (6.2), 
we simulated m values ranging from 0.00375 to 0.25 in lattice units and searched numerically for a simple 
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scaling form, Eq. (6.2). We found that very accurate measurements of E were possible. Vacuum tunnelling 
was always suppressed in these measurements because of the explicit symmetry breaking in the action. Our 
computer results are given in Tables X and XI. 

In Fig. 15 wc plot the data for a sequence of lattices 8 3 , 16 3 and 24 3 . In fact, In So is plotted against 
lnm for 1/g 2 — 1.00. We see finite size effects at the smallest m values as expected - as m decreases larger 
and larger lattices are needed to measure S 's nonzero value. For m ranging from 0.0625 to 0.01625 we see 
the 8 3 and 16 3 data approaching the 24 3 data from above. A linear fit over this range of m gives 8 = 2.0 in 
almost perfect agreement with the large Nf predictions. Note that fermion masses less than 0.01625 were 
not used in the fit because the 16 3 and 24 3 lattices do not give consistent values for S , indicating that finite 
size effects arc still important here. Presumably, if one simulated 32 3 lattices, a wider "scaling window" 
would be seen. Note that error bars on the S measurements are smaller than the symbols in the figure. 

In Fig. 16 we show identical plots for simulation runs at 1/g 2 = 0.975 which could also be the bulk 
critical point. The value obtained for 8 for the fit here is slightly larger, 2.21, than that discussed above. We 
again see that without better knowledge of the bulk critical coupling our uncertainty in the critical index 
8 will be about 10%. Uncertainty in the position of critical points usually limits the precision of finite size 
scaling studies in other statistical mechanics models as well. 

Finally, we measured 8 for the finite temperature rather than the bulk critical point. This is particularly 
interesting because the large Nf prediction for 8 associated with the finite temperature transition is 3 [17], 
indicative of a Gaussian model rather than the nontrivial scaling theory at the bulk critical point. To 
determine this critical index we used the largest lattice studied in Sec. 4, 12 x 36 2 and added a bare mass 
term to the action. In this case the critical coupling for the finite temperature transition on the 12 x 36 2 
lattice was determined accurately in Sec. 4, l/g 2 c (N T = 12) = 0.880(5), so a relatively high precision 
measurement of its 8 is possible. The computer data are recorded in Tabic XII. In Fig. 17 we show lnS 
vs. lnm for m values ranging from 0.01 to 0.000625. Linear fits to the data produce the result 8 = 3.0(1) in 
excellent agreement with our analytic expectations. 

7. The Scalar Correlation Function 

In this section we report on measurements of the scalar two-point function < a(0)<r(x) >, and our 
attempts to understand them using the predictions of the 1/Nf expansion described in Sec. 2. Since we 
have been interested for the most part in the bulk properties of the model, we did not perform any special 
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simulations on lattices which are longer in one direction, which is normal practice in numerical QCD, where 
the requirement is to examine the asymptotic behavior of hadronic correlators in order to extract hadron 
masses. Further studies may benefit from this approach: however, as we shall see, the relation between the 
asymptotic decay of the correlator and the particle mass is by no means so straightforward in this model, and 
indeed, in the symmetric phase we are only able to extract a width by fitting to the sub-asymptotic behavior! 
Instead, we performed measurements on 20 3 lattices at the values of the coupling 1/g 2 used elsewhere. 

As usual, in order to simplify the fit and maximize statistics we projected on to the zero spatial momentum 
sector by summing over timeslices, and define the correlator: 



where y labels the V s sites existing on each of N T timeslices. Of course, the quantity of interest is really the 
connected correlator 



In our work we assume the expectation value < a > is translationally invariant, and extract < er > 2 as a 
fitted parameter in fits of C (x) . This is a much more stable procedure than subtracting the measured value 
Sq from the data before fitting, because the vacuum and the other quantities entering the fit are highly 
correlated. Of course, we expect agreement between the two quantities. 

In order to parameterize the decay of C c (x), we use the leading order expressions for D a given in Sec. 2. 
First we discuss fits in the broken phase, where the scalar is expected to be a genuine massive bound state. 
From (2.7) we have 



where here M denotes the scalar mass. Immediately we see the difference between this model and QCD: the 
presence of the hypergeometric function in the denominator introduces branch cuts in the complex fc-planc 
along the imaginary fc-axis in the ranges [iM,ioo) and (— ioo, —iM]. Instead of a discrete series of poles, 
corresponding to a set of bound states in the channel in question, which is the case in a confining theory 
such as QCD, we have a complicated spatial dependence due to a strongly-interacting fermion - anti-fermion 
continuum for k > M a = 2S (to leading order in l/Nf) [9]. Performing the contour integral around the 
branch cut in the upper half-plane, we find 




(7.1) 




(7.2) 




(7.3) 




(7.4) 
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where the second form is suitable for numerical quadrature. Due to the contribution of the fermion - anti- 
fermion continuum, P(x;M) decays faster than e~ Mx , and the corrections to the latter form cannot be 
expressed as a discrete sum of exponentials over higher states. This means that we no longer have methods 
such as a search for a plateau in the effective mass plot to tell us when the form (7.4) might become valid: 
we must simply trust its validity over the whole range of x, or at least that range exceeding the correlation 
length, where we might expect the continuum approximation to be accurate. Note that in ref. [31] the scalar 
correlator in the two-dimensional Gross-Neveu model, which is asymptotically free, was analyzed in this way, 
and it was found that the inclusion of a free fermion continuum in the fit did improve the results. In our 
case, however, the fermions are strongly interacting, so we are forced to try (7.4). 

We have fitted C(x) for l/g 2 ranging from 0.7 to 0.975, using the three-parameter form a[P(x; M) + 
P(L — x; M)\ + £q, where L = 20 is the length of the lattice. The results are presented in Table XIII and 
plotted in Fig. 18 (error bars are not shown for clarity). We used a least squares fit with a correlated 
X 2 function [33] to allow for correlation between timeslices. Errors were estimated using standard binning 
procedures. For comparison we also tried a standard single pole fit of the form a[e~ Mx + e~ M ( L ~ x *>] + Eg. 
In both cases we fitted for x ranging between 2 and 9; this produced % 2 values per degree of freedom < 1 in 
most cases (the maximum x 2 per d.o.f. found was 1.64). 

We can make the following observations. The values of M obtained by the branch cut form P(x; M) 
lie systematically below those obtained from the single pole form, as expected: however the fluctuations in 
the two are clearly correlated. We see that M decreases as g 2 — > g 2 + ~ 0.96 on the 20 3 lattice, but does 
not reach zero as required by the hypothesis that it represents an inverse correlation length. The points at 
l/g 2 = 0.95, 0.975 clearly lie outside this trend, probably because they are in or very close to the symmetric 
phase. The scatter of the M values about the downward trend is accounted for by the errors given in Table 
XIII. We also plot the fitted values of 2Eg (which coincide within error bars for both fits), which at leading 
order in 1/Nf corresponds to the two fermion threshold (note that even for a correlation length of 10, the 
0(1/Nf) corrections predicted by Eq.(2.15) amount to just 5% for Nf = 12). The line drawn through these 
points is that obtained from a fit to the data of Table VI, so we see the fitted values of Eg coincide with 
good accuracy with the bulk averages, which is a consistency check of our fitting method. 

If the scalar field interpolates a true bound state, then M < 2E . We expect the difference, representing 
the binding energy, to be 0(1/Nf). We see that our fits for M using P(x;M) obey this criterion only for 
l/g 2 < 0.875; moreover there is no clear region of the plot where the ratio M : 2S is constant, so no estimate 
of M in physical units can be made. Given the systematic uncertainties in the fitting procedure highlighted 
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by the difference between the two fitted forms, we interpret this as an indication that the leading order 
form (7.3,4) is still not adequate to fit the measured correlator. Two further avenues suggest themselves. 
First, one could calculate the leading order scalar propagator in an explicit lattice regularization: this would 
probably need to be done numerically, but might, for instance, modify parameters such as the tt 2 in the 
denominator of (7.4) which at present are put into the fits "by hand" . Secondly, one could bite the bullet 
and attempt to calculate D (T (k 2 ) to 0(1/Nf) in the broken phase, in order to be able to parameterize the 
separation between the bound state pole and the two-fermion threshold. Our attempts to do this in a naive 
way did not yield useful results. 

Next we consider the scalar correlator in the symmetric phase. As discussed in Sec. 2, in this case 
the scalar field does not interpolate a stable massive particle state, but rather an unstable resonance. The 
resonance width fi serves as an inverse correlation length, and can be extracted from the data as follows. 
We start from the definition (7.1) of C(x) and the leading order form for D a (k 2 ) given in (2.26): 

C(x) oc / = Re (e-^EA-ifix)) = g(fix). (7.5) 
Jo k + fi 

Here E\(z) is the exponential integral function defined by [34]: 

f°° e _t 

Eh.(z) = J —dt. (7.6) 



fW = 732 i + fe))' ( 7 - 7 ) 



For fix — > oo we have 

In this phase the scalar correlator decays asymptotically with a power law: the corrections to this power 
law at small fix enable us to extract the physical scale fi. Because of the power-law form, it is not sufficient 
to sum over just forward and backward propagating signals to get an accurate fit: one must sum over the 
signals from all "images" corresponding to propagation an arbitrary number of times around the periodic 
lattice. Using the asymptotic form (7.7) and the Euler-McLaurin formula, we can approximate this sum: 



C(x) cx Q(x; fi) = ^2 g(fi(x + nL)) 

n— — oo 

= ± fl(M(a + ^) ) + ^( ] ^ + _!_) 

n— — N 



(7.8) 



2fi 2 \(NL + x) 2 (NL-x) 2 J 6fi 2 \(NL + x) 3 (NL-x 

In Table XIV and Fig. 18 we show the results of the two-parameter fit C(x) — aQ(x;fi) to results 
obtained for 1/g 2 between 1.1 and 0.925. Once again, we fitted to x values 2-9, and the % 2 per d.o.f. was 
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in most cases less than one, with a maximum value 1.58. Fits of the standard simple pole form in this region 
were of worse quality and yielded higher values of \i: moreover if we constrained the vacuum S = 0, then 
these fits became disasterously poor. In most cases we found the series truncation in (7.8) was adequate for 
N = 5. Only for the very smallest \i at 1/g 2 — 0.925 did we find some weak dependence of the fit on N - 
here we quote the N = 40 result. 

From Fig. 18 we see that in the symmetric phase the fitted parameter \x decreases fairly smoothly to 
zero as g 2 — > g 2 _, and has every appearance of a genuine inverse correlation length. The data are not yet 
of sufficient quality to make a meaningful extraction of the exponent v, however. It is also not clear how 
reliable the points at 1/g 2 = 0.925 and 0.95 are - perhaps a study on different sized lattices might be useful 
here. Overall, however, given the limitations of our statistics and lattice sizes, we find the results of the fit 
quite satisfactory. This is probably a sign that the leading order form of the propagator in the symmetric 
phase (2.26) is not grossly changed by 1/Nf corrections: this is of course verified by (2.35), where we find 
the only change to D a (k 2 ) is an 0(1/Nf) correction to the power of k. 

To sum up, we have identified a correlation length associated with the auxiliary field a in both phases 
of the model, which increases as 1/g 2 — > l/g 2 ±- Leading order 1/Nf calculations predict that the inverse 
correlation length should decrease linearly in this limit in either phase, ie. v = 1+0(1/Nf). In the symmetric 
phase our accuracy appears limited principally by statistics and lattice volume; the leading order prediction 
is qualitatively verified, but we are unable to fit directly for the exponent v. In the broken phase, the 
correlation length we extract from our fits does not diverge, and we suspect an improved understanding of 
the 1/Nf or O(a) corrections in this phase will be required before more progress can be made. 

8. Summary 

In this article we have attempted to reconcile two alternative views of the strongly-coupled critical point 
which is generic to the model four-fermi Lagrangian (1.1). In the language of statistical mechanics, we have 
characterized the fixed point by calculating its critical exponents, which in turn define its universality class. 
The indices satisfy consistency conditions known as hyperscaling relations, which supports the idea of a single 
divergent correlation length in the critical region. The critical correlation functions display non-canonical 
scaling, that is, the operators corresponding to physical observables acquire anomalous dimensions. 

In particle physics language, we have seen that for some sufficiently large coupling, the most stable 
vacuum is one in which chiral symmetry is spontaneously broken, leading to a dynamically-generated fermion 

38 



mass. An expansion in powers of l/Nf about the critical point is exactly renormalizable for a continuum 
of dimensionality 2 < d < 4 - this behavior is unlike standard weak coupling perturbation expansions 
(WCPE) about the Gaussian fixed point, which are exactly renormalizable only for a single critical value 
of d. Interactions in the model occur via exchange of scalar degrees of freedom. The composite nature of 
these particles ensures that the physical interaction strength remains non- vanishing as the UV cutoff is sent 
to infinity - the theory is non-trivial. 

We have argued, in general terms in section 1 and in more detail in section 2, that these descriptions are 
equivalent, and that rcnormalizability, hyperscaling and non-triviality are inextricably linked in this model. 
The composite nature of the a field is signalled by its large anomalous dimension. The main results of the 
explicit calculations of 0(1/Nf) corrections in section 2 are as follows: firstly, at this order and beyond, 
logarithmic divergences appear, which necessitate vertex, wavefunction, and mass renormalizations. The 
set of primitively divergent diagrams is such that the renormalization constants are overdetermined, which 
means that renormalizability depends on a non-trivial cancellation between different one-particle irreducible 
(1PI) vertex functions. This cancellation can be viewed as a constraint on the critical indices equivalent to 
hyperscaling, but ultimately results from the existence of Ward identities arising from the chiral symmetry 
of the model. Secondly, we have seen that the critical exponents get explicit l/Nf corrections, implying 
that theories with different Nf belong to different universality classes. As we calculate to higher order in 
l/Nf, we improve our understanding of the fixed point theory - this is in contrast to WCPE, where the fixed 
point is known to be Gaussian from the start, and perturbative corrections serve to parameterize trajectories 
in or near the critical hypersurface, that is, to relate massless theories at differing momentum scales. In 
the latter case we say the WCPE interaction is marginal, whereas in the four-fermi model the interaction 
(■tpip) 2 is relevant, since any departure from the critical coupling strength g 2 ^ g 2 results in the theory 
becoming massive, with a finite correlation length, either in the conventional sense (t > 0), or by generating 
a dimensionful width for the a resonance (t < 0). So far we have no reason to suspect the critical hypersurface 
has a dimension greater than zero [22]. In simple terms, the l/Nf expansion can never yield a "running 
coupling constant", since Nf is a basic parameter of the theory and is not itself renormalized. Thirdly, 
as a result of the consistency between the different 1PI amplitudes, or if preferred the chiral symmetry of 
the model, we see that in the deep Euclidean limit k 2 — > oo the four-fermi scattering amplitude has no 
dependence on l/Nf, apart from a trivial overall factor, and assumes a universal form Ad/N fk d ~ 2 . This 
is the universal interaction which characterizes the UV fixed point. It is interesting that despite the fact 
that we have used the l/Nf expansion to calculate radiative corrections, the l/Nf corrections manifest 
themselves in the IR properties of the theory, such as the bulk critical exponents, and presumably the bound 
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state spectrum, about which we have little knowledge at present. Finally, we have shown that the picture 
we have developed is readily extended to models with a continuous chiral symmetry, so that there appear to 
be no problems incorporating massless Goldstone excitations. 

On the numerical side, our main result is that the leading order predictions of the 1 /Nf expansion have 
been verified in a non-perturbative simulation, suggesting that the 1/Nf expansion is accurate at least for 
Nf > 12. The critical exponents (3, 7, 5 and v have all been directly estimated and found to be distinctly 
non-Gaussian - this is strong evidence for the presence of anomalous dimensions. Moreover, simulations at 
non-zero temperature support a crossover from the leading order result 5 = d — 1 to the Gaussian mean 
field value (5 = 3. We have not been able to measure explicit 1/Nf corrections in this study, but are able 
to estimate the extra computational effort required. We have also laid the groundwork for simulations at 
smaller Nf, where quantum fluctuations will be larger, and the 1/Nf corrections easier to detect. It is, 
of course, possible that for some sufficiently small Nf the picture we have developed breaks down and the 
1/Nf expansion loses its applicability - one suggestion has been that the transition becomes first order 
[35]. Another direction, which we have explored in a separate paper [32], is to simulate the model at non- 
vanishing chemical potential. Unlike the case for lattice gauge theories, the action (3.1) remains real under 
this modification, which means that we are able to study a relativistic theory of interacting fermions at 
non-zero density on far larger systems than are otherwise possible. As reported in [32], these simulations 
have enabled an independent estimate of the exponent v. 

Progress on understanding the scalar two-point function has been slower, highlighting the difficulties in 
understanding the spectrum of a theory which is both strongly coupled and non-confining. The main novelty 
and success here was in the symmetric phase, where a "mass" (inverse correlation length) was extracted from 
a fit to the pre-asymptotic behavior of the correlator. Further progress in this direction probably depends 
upon a more complete understanding of 0(1 /Nf) corrections in the broken phase, which is technically much 
tougher than the calculation presented here. Of course, in the long run one would wish for a complete 
understanding of the bound-state spectrum. 

We hope that many of the aspects we have learned about this model will be of relevance to strongly 
coupled quantum electrodynamics in four dimensions, which has been the focus of much recent work by 
ourselves and others. Both models exhibit a continuous phase transition from a "massless" phase to one 
in which mass is dynamically generated - and both models appear to scale with non-Gaussian critical 
exponents [36]. The great hope is, of course, that QED4 will turn out to be the first known example of a 
non-trivial non-asymptotically free theory in four dimensions. As yet, we have only analytic arguments based 
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on self-consistent approximations, and numerical simulations, which, however powerful, are always in need 
of some kind of interpretation. In the models discussed here we have the additional benefit of a systematic 
approximation, which has enabled us to make precise calculations and gain physical insight. Perhaps the 
ultimate result of this work will be the development of a concise vocabulary to describe strongly- coupled UV 
fixed points. 



Here we give some technical details of how to perform various momentum integrals in d-space. First, we 
define our Dirac matrices for all values of d to be 4 x 4, traceless, and hermitian, so that the trace always 
yields a factor of 4. Alternative definitions are possible [21], but the results are qualitatively unchanged, at 
least to 0(1 /N f ). 

We define the integral over Euclidean d-momcntum p as follows: 



Although the formulae are superficially similar, this is not dimensional regularisation; the cutoff A appears 
explicitly in (A.l) and is used to regulate divergent integrals. The following formula, which holds for 
d G (2,4), serves to evaluate loop integrals involving only fermi propagators: 



Finally, we give a brief outline of how to evaluate the two-loop integral (2.29). With the momentum 
routing specified the integral over p is fairly straightforward: after taking the trace it reduces to several 
components which can each be brought to the form (A. 2) using a suitable Feynman parameterization and 
momentum shift. If only one Feynman parameter is required in a particular sub- integral, it is then simple 
to extract divergent terms of the power-law form (2.30). If there are two or more Feynman parameters, then 
there will be an intermediate integral with a denominator, which is a function of q and external momentum, 



Appendix 




(A.l) 
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raised to a non-integer power s. The following formula is needed: 



/ 



(a + bx + cx 2 ) s \ 2c 2c 2 / V A / V 2c 

bC\ (4c) 8 - 1 



2(1 -s)c A + 4c 2 (.t + 



2c; 



(A3) 



CA /4c V / 6 \ „ / , , 4c 2 



4c 2 V A/ V 2c 



^(-i.J;i;-x( a! + =) a ) 



where A = 4ac— b 2 and F is the hypergeometric function. This result is exact. However, it turns out that the 
momentum-dependent logarithmic divergences (2.31) naturally arise from the lower limit of the Feynman 
parameter integration, ie. from setting x = above. Also, in practice, the coefficients b, c are 0(q 2 ), where q 
is the remaining loop momentum, but a is 0(q°). These two considerations make the following simplifications 
possible: 

J Q (a + bx + cx 2 Y = ~ A{s - l)^- 1 + ° [y) 5 (A4fl) 
^ - 1 (3 - 2s) + 0(\\ (AAb) 



(a + bx + cx 2 ) s A(s-l)a s " 2 A(s - 2)(s - l)a s - 2 V 9' 



i 



2<3 -" )b + *Zj*+o(l). (A. 4c) 



l (a + bx + cx 2 ) s A 2 (s - 2)(s - l)^- 3 A 2 (s - 3)(s - 2)a s - 3 

With judicious choice of Feynman parameterization these formula; suffice to evaluate (2.29). For example, 
suppose after taking the trace, introducing Feynman parameters x and y, and shifting the momentum p, we 
are left with the following p sub-integral: 

I(q,k) = 8 [ dx [ X dy [ — - xk 2 + yk.q 

1 ' Jo Jo J P [p 2 +x{l-x)k 2 +y{l-y)q 2 -2xyk.qY V ' 

We can do the integration over p using (A. 2): 

rt is , r ( 3 -!) (\ f 1 '* , xk 2 +yk.q . . 

I(q,k)=4— / dx / dy y H — r . (A6) 

(47r)2 Jo Jo [x(1 - x)k 2 + (q 2 - 2xk.q)y - q 2 y 2 ] * 

Now use (A. 4) to do the integral over y; only the lower limit is relevant: 

/(g,fc)^4 r(3 "| ) C dxxi-'il-xf^^^i-, \ +o(\)). (A.7) 

V ^ (4tt)* Jo V ' {2-±)\q 2 + Ax{l-x)k 2 WJJ 

The remaining integral over x, and the subsequent integral over q are now straightforward, the final contri- 
bution to (2.29) being 

Z a g 2 r(2-l)B(l,l-l)A d A 
-~ 8 *T (4^rpT(|) (M 
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Table I 



Table of critical exponents and critical couplings 1/G 2 = 1/gl x (47r)2r(|)(d - 2)/8A d 2 for the three 
variants of the Gross-Neveu model considered in section 2. The constant C d is defined in Eq.(2.13). 
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4-d + 



(4-d) C d 
d N f 



1 



27V/ 



(^-2) 
N f 



(2d - 5) 



2 AT/ 



4G 



Table II 



Vacuum expectation value £o(iV/ = oo) as predicted by the leading order lattice gap equation, together with 
measured deviations A r /(S (oo) — Eo(iV/)) vs. coupling 1/g 2 , for Nf = 6, 12,24. 



1/g 2 EoM 24(S (oo)-S (24)) 12(E„(oo) - S (12)) 6(S (oo) - E (6)) 



.5 


.8351 


.199(14) 


.227(12) 


.267(10) 


.6 


.6401 


.310(19) 


.319(12) 


.382(11) 


.7 


.4711 


.408(22) 


.402(13) 


.517(15) 


.8 


.3165 


.542(34) 


.550(24) 


.634(30) 
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Table III 



Vacuum expectation value S an d its susceptibility \ vs. coupling 1/g 2 on a 8 3 lattice. The number of 
trajectories of the hybrid Monte Carlo algorithm at each value of 1/g 2 are listed in column 4. 



1/g 2 Eo X Trajectories 



.70 


.4397(12) 


.265(8) 


5,000 


.725 


.3974(14) 


.308(12) 


5,000 


.75 


.3544(16) 


.400(14) 


5,000 


.775 


.3114(20) 


.534(25) 


5,000 


.80 


.2693(20) 


.721(31) 


7,500 


.825 


.2136(21) 


1.485(35) 


7,500 


.85 


.1579(19) 


2.07(14) 


10,000 


.875 


.1152(19) 


1.75(18) 


10,000 


.90 




1.38(18) 


10,000 


.925 




.916(44) 


10,000 


.95 




.633(41) 


7,500 


.975 




.489(32) 


7,500 


1.00 




.363(21) 


5,000 


1.025 




.270(9) 


5,000 


1.05 




.215(7) 


5,000 
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Table IV 



Vacuum expectation value So and its susceptibility \ vs - coupling 1/g 2 on a 12 3 lattice. 



1/g 2 Eo X Trajectories 



.70 


.4319(10) 


.511(11) 


5,000 


.725 


.3905(12) 


.621(13) 


5,000 


.75 


.3482(12) 


.695(14) 


5,000 


.775 


.3091(13) 


.876(25) 


5,000 


.80 


.2655(15) 


.974(33) 


5,000 


.825 


.2255(15) 


1.06(7) 


7,000 


.85 


.1786(17) 


2.34(17) 


7,000 


.875 


.1262(21) 


3.28(75) 


10,000 


.90 


.0853(25) 


3.38(81) 


10,000 


.925 




6.38(2.00) 


10,000 


.95 




3.75(1.00) 


10,000 


.975 




2.18(22) 


5,000 


1.00 




1.87(21) 


5,000 


1.025 




1.53(16) 


5,000 


1.05 




1.22(4) 


5,000 
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Table V 



Vacuum expectation value So and its susceptibility \ vs - coupling 1/g 2 on a 16 3 lattice. 



1/g 2 Eo X Trajectories 



.70 


.4323(15) 


.546(24) 


1,000 


.725 


.3888(15) 


.556(25) 


1,000 


.75 


.3462(17) 


.741(27) 


1,000 


.775 


.3065(18) 


.826(45) 


1,000 


.80 


.2625(21) 


.988(60) 


1,000 


.825 


.2219(22) 


1.19(8) 


1,000 


.85 


.1747(23) 


1.85(30) 


2,000 


.875 


.1306(25) 


3.05(35) 


2,000 


.90 


.0885(35) 


4.12(1.12) 


2,000 


.925 




10.6(4.0) 


2,000 


.95 




4.97(1.60) 


2,000 


.975 




2.94(50) 


1,000 


1.00 




3.22(50) 


1,000 


1.025 




1.22(14) 


1,000 


1.05 




.97(8) 


1,000 
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Table VI 



Vacuum expectation value So and its susceptibility x vs. coupling 1/g 2 on a 20 3 lattice. 



1/g 2 Eo X Trajectories 



.70 


.4317(2) 


.611(15) 


5,000 


.725 


.3883(2) 


.633(15) 


5,000 


.75 


.3460(2) 


.746(25) 


5,000 


.775 


.3045(3) 


.837(51) 


5,000 


.80 


.2616(3) 


.948(55) 


7,500 


.825 


.2189(3) 


1.21(8) 


7,500 


.85 


.1800(4) 


1.40(9) 


7,500 


.875 


.1376(5) 


2.24(15) 


10,000 


.90 


.0873(15) 


4.07(1.15) 


10,000 


.925 


.0406(20) 


3.86(1.16) 


10,000 


.95 




7.49(2.50) 


7,500 


.975 




4.19(1.20) 


7,500 


1.00 




2.28(32) 


7,500 


1.025 




1.66(15) 


5,000 


1.05 




1.35(10) 


5,000 


1.075 




1.05(7) 


5,000 


1.10 




.770(2) 


5,000 
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Table VII 



Estimates of the bulk critical point on L 3 lattices with L = 8, 12, 16, 20. 



L i/g 2 c (L) 

8 .867(8) 

12 .925(20) 

16 .930(20) 

20 .950(20) 



Table VIII 

Estimates of the finite temperature critical points on N T x N 2 lattices for N T = 2, 4, 6, 8, 10 and 12. 







2 


.47(1) 


4 


.695(5) 


6 


.785(5) 


8 


.833(5) 


10 


.865(5) 


12 


.880(20) 
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Table IX 



Zero temperature vacuum expectation values of S (from Fig. 9) at the couplings of the finite temperature 
critical points recorded in Table VII, and T c estimates for N T ranging from 2 through 12 in units of So- 





So 


T c /E 


2 






4 


.44(1) 


.57(1) 


6 


.29(1) 


.57(2) 


8 


.21(1) 


.60(3) 


10 


.157(10) 


.64(4) 


12 


.125(10) 


.67(5) 
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Table X 



Vacuum expectation value S vs. fermion bare mass m at coupling 1/g 2 — 1.00 on 8 3 , 16 3 and 24 3 lattices. 
The number of trajectories used for each measurement are given in the last column. 



TO 




S (16' 3 ) 


S (24 J ) 


Trajector 


.005 


.0496(12) 


.0614(10) 


.0594(10) 


10,000 


.008 


.0703(13) 


.0783(9) 


.0814(10) 


10,000 


.0125 


.0931(5) 


.106(1) 


.106(1) 


10,000 


.01625 


.118(2) 


.122(1) 


.123(1) 


10,000 


.020 


.138(1) 


.138(1) 


.137(1) 


10,000 


.025 


.155(1) 


.154(1) 


.153(1) 


10,000 


.0375 


.190(1) 


.187(1) 


.187(1) 


10,000 


.050 


.218(1) 


.213(1) 


.213(1) 


5,000 


.0625 


.239(1) 


.234(1) 


.235(1) 


5,000 


.075 


.257(1) 


.253(1) 


.252(1) 


5,000 


.10 


.287(1) 


.282(1) 


.282(1) 


5,000 


.15 


.327(1) 






5,000 


.20 


.355(1) 






5,000 


.25 


.377(1) 






5,000 
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Table XI 



Vacuum expectation value S vs. fermion bare mass m at coupling 1/g 2 = 0.975 on 8 3 , 16 3 and 24 3 lattices. 



m 



So(8 3 ) 



S (16 3 ) 



So(24 : 



Art r 

.005 


.0647(10) 


DTD E / C \ 

.0705(5) 


.0767(3) 


10,000 


.008 


.0943(8) 


.0968(3) 


.0980(3) 


10,000 


.0125 


.118(2) 


.123(1) 


.122(1) 


10,000 


.01625 


.134(1) 


.135(1) 


.129(1) 


10,000 


.020 


.150(1) 


.152(1) 


.152(1) 


10,000 


.025 


.170(1) 


.170(1) 


.169(1) 


10,000 


.0375 


.207(1) 


.203(1) 


.202(1) 


10,000 


.050 


.233(1) 


.229(1) 


.228(1) 


5,000 


.0625 


.254(1) 


.250(1) 


.249(1) 


5,000 


.075 


.272(1) 


.267(1) 


.267(1) 


5,000 


.10 


.300(1) 


.297(1) 


.297(1) 


5,000 



Trajectories 
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Table XII 



Vacuum expectation value £ vs. fermion bare mass m on a 12 x 36 2 lattice at coupling 1/g 2 — 0.8775. The 
number of trajectories used for each measurement are given in the last column. 



m £ Trajectories 



.01 


.1826(8) 


1,000 


.005 


.1471(4) 


1,000 


.0025 


.1151(2) 


1,000 


.0018 


.1084(2) 


1,000 


.00125 


.0903(10) 


1,200 


.0009 


.0805(10) 


1,200 
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Table XIII 

Results of fits to the scalar correlator in the broken phase, using both the branch cut form P(x; M) and the 
usual simple pole form. 







branch cut 






simple pole 




v y 


M 





A. 


M 


S 2 , 


A 


.7 


.708(98) 


.18606(1) 


3.3 


.812(94) 


.18606(1) 


2.9 


.725 


.750(95) 


.15069(1) 


4.2 


.836(87) 


.15069(1) 


4.3 


.75 


.494(69) 


.11958(1) 


3.4 


.607(64) 


.11958(1) 


3.7 


.775 


.388(58) 


.092382(18) 


3.3 


.513(51) 


.092391(15) 


5.1 


.8 


.433(58) 


.068532(17) 


6.2 


.553(51) 


.068539(15) 


5.9 


.825 


.348(48) 


.047656(23) 


8.2 


.478(40) 


.047667(17) 


6.7 


.85 


.320(57) 


.032503(30) 


2.1 


.450(50) 


.032519(22) 


1.9 


.875 


.290(36) 


.018956(25) 


3.1 


.433(27) 


.018982(16) 


5.1 


.9 


.281(42) 


.007620(32) 


4.7 


.428(29) 


.007653(19) 


2.5 


.925 


.247(55) 


.001906(46) 


3.1 


.387(39) 


.001941(23) 


2.4 


.95 


.354(57) 


.000812(24) 


2.0 


.475(51) 


.000821(19) 


2.9 


.975 


.312(61) 


.000383(29) 


3.1 


.459(51) 


.000407(19) 


2.5 
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Table XIV 

Results of fitting the scalar correlator in the symmetric phase to the form Q(x; ji). 



1 /«2 


/i 


X 


1.1 


.437(55) 


5.8 


1.075 


.283(33) 


3.1 


1.05 


.159(14) 


9.5 


1.025 


.156(13) 


4.6 


1.0 


.146(12) 


4.0 


0.975 


.075(5) 


4.7 


0.95 


.041(2) 


8.3 


0.925 


.019(1) 


4.8 
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Figure Captions 

Figure 1: 0(1 /Nf) corrections to the fermion self-energy and the fermion-scalar vertex. 
Figure 2: 0(1 /Nf) correction to the gap equation. 
Figure 3: 0(1 /Nf) corrections to the scalar propagator. 

Figure 4: Schematic diagram showing 1/Nf corrections to the four-fermion scattering amplitude. 

Figure 5: Plot of S vs. 1/g 2 on a 12 3 lattice, for Nf = 24 (squares), Nf = 12 (circles), and Nf = 6 
(triangles). The solid line is the leading order solution to the lattice gap equation (3.13). Errors are smaller 
than the size of the symbols. 

Figure 6: Plot of the data recorded in Table III. the reciprocal of the susceptibility, is plotted multiplied 
by 10" 1 . 

Figure 7: Plot of the data recorded in Table IV. x _1 is plotted multiplied by 10 _1 . 

Figure 8: Plot of the data recorded in Table V. x _1 is plotted multiplied by 10 _1 . 

Figure 9: Plot of the data recorded in Table VI. x _1 is plotted multiplied by 10 _1 . 

Figure 10: Plot of m£ vs. ln(l/# 2 - i/g*) f rom Table VI, 1/g 2 = 0.955. 

Figure 11: Plot of the bulk critical point's dependence on lattize size, l/g 2 (L) vs. 1/L. 

Figure 12: Histograms of £ measurements on a 10 x 30 2 lattice for couplings 1/g 2 = 0.870, 0.865 and 0.860. 

Figure 13: Plot of the finite temperature critical point's dependence on the temporal extent N T of the 
asymmetric N T x N 2 lattice. 

Figure 14: Plot of ln(l/^ - l/g 2 c (N T )) vs. lnN T for three estimates (0.995, 0.976, 0.950) of the infinite 
volume critical point. 

Figure 15: Plot of ln£ vs. In to for 8 3 , 16 3 and 24 3 lattices at 1/g 2 = 1.00, an estimate of the bulk critical 
point. 

Figure 16: Plot of lnS vs. In to for 8 3 , 16 3 and 24 3 lattices at 1/g 2 = 0.975, an estimate of the bulk 
critical point. 

Figure 17: Plot of InS vs. In to on a 12 x 36 2 lattice at 1/g 2 = 0.8775, the finite temperature critical point 
for this lattice size. 

Figure 18 Plot of fits of the inverse correlation length vs. 1/g 2 for a 20 3 lattice. The pluses show the scalar 
mass M in the broken phase obtained using a fit to a branch cut (7.4), the stars show M values obtained 
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using a simple pole fit, and the circles show the 2 fermion threshold 2S obtained using the branch cut fit. 
In the symmetric phase the crosses show the scalar width \i fitted using the form (7.8). 



GO 



